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ABSTRACT 



A transient three-dimensional heat conduction code was devel- 
oped using finite differences. A stability restriction on the time step 
was avoided using a technique proposed by Brian. Computations from 
the code were validated using both the explicit technique and an 
available closed-form solution for small times. The maximum error was 
found to be within 0.019 percent for an llxllxll grid and time 
step of 17.117 seconds. The total CPU time to carry out the computa- 
tions up to 3,600 seconds using Brian’s technique was six times that 
required for the explicit technique with the same time step of 17.117 
seconds. However, as the time step was increased without altering the 
geometry, the CPU time using Brian’s technique decreased and was 
less than that used in the explicit technique for time steps larger than 
110 seconds. The validated code was also used in the analysis of the 
transient thermal response of a component on an orbiting spacecraft. 
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THESIS DISCLAIMER 

The reader is cautioned that computer programs developed in 
this research may not have been exercised for all cases of interest. 
While every effort has been made, within the time available, to ensure 
that the programs are free of computational and logic errors, they 
cannot be considered validated. Any application of these programs 
without additional verification is at the risk of the user. 



iv 



TABLE OF CONTENTS 



I. INTRODUCTION 1 

A STATEMENT OF PROBLEM 1 

B OBJECTIVES 2 

II. DEVELOPMENT OF MODEL 3 

A MODEL CONDITIONS 3 

B EXPLICIT TECHNIQUE 4 

1 . Derivation of Nodal Temperature Equations 4 

2. Restrictions on Time Step in the Explicit Technique 7 

3. Numerical Calculations 9 

C BRIAN’S TECHNIQUE 9 

1. Derivation of Nodal Temperature Equations 1 1 

2. Solution of the Nodal Temperature Equations 1 7 

3. Numerical Calculations 1 8 

III. VALIDATION 20 

A SELECTED GEOMETRY AND CONDITIONS 20 

B NUMERICAL COMPUTATIONS 22 

C VALIDATION COMPUTATIONS 24 

IV. APPLICATIONS 35 

A GENERAL 35 



v 



R ORBITING SATELLITE COMPONENT 35 

1. Geometry and Initial Conditions 35 

2. Numerical Calculations 4 1 

V. RESULTS 45 

VI. CONCLUSIONS 46 

VII. RECOMMENDATIONS 47 

APPENDIX A EXPLICIT TECHNIQUE NODE EQUATIONS 48 

APPENDIX B BRIAN’S TECHNIQUE NODE EQUATIONS 57 

APPENDIX C PROGRAM EXPLICIT 8 5 

APPENDIX D PROGRAM BRIAN 1 0 7 

APPENDIX E PROGRAM VALID 1 8 7 

LIST OF REFERENCES 191 

INITIAL DISTRIBUTION LIST 192 



vi 



LIST OF TABLES 



1 Parameters Used in Validation Runs 2 1 

2 Parameters Used in Orbiting Satellite Application 40 



vii 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 1 

12 

13 

14 



LIST OF FIGURES 



Left Surface of Rectangular Block Used in the Model 5 

Flow Chart for Program EXPLICIT 1 0 

Flow Chart for Program BRIAN 19 

Rectangular Model Showing Boundary Conditions Used in 
Validation Runs 2 1 

Flow Chart for Program VALID 2 3 

Temperature Difference Versus Time for Time Increment of 

17.117 Seconds with Distance of 1 Meter Into Block 25 

Temperature Difference Versus Time for Time Increment of 

4.279 Seconds with Distance of 1 Meter Into Block 26 

Temperature Difference Versus Time for Time Increment of 

17.117 Seconds with Distance of 0.4 Meters Into Block 27 

Temperature Difference Versus Time for Time Increment of 

4.279 Seconds with Distance of 0.4 Meters Into Block 28 

Temperature Difference Versus Time for Time Increment of 

17.117 Seconds on Left Surface of Block 2 9 

Temperature Difference Versus Time for Time Increment of 

4.279 Seconds on Left Surface of Block 30 



Temperature Difference Versus Time for Time Increments of 

17.117 and 120 Seconds with Distance of 1 Meter Into Block.... 31 



Temperature Difference Versus Time for Time Increments 

of 17.117 and 120 Seconds with Distance of 0.4 Meters 

Into Block. 32 

Temperature Difference Versus Time for Time Increments 
of 17.117 and 120 Seconds on Left Surface of Block 33 



viii 



15 Component on a Satellite in Circular Orbit 37 

16 Geometry Used in Determining Geometric Factor for Earth 

Thermal Radiation 38 

1 7 Geometry Used in Determining Geometric Factor for Solar 

Reflected Radiation 39 

1 8 Temperature Variation of the Top Surface Center Node 42 

1 9 Temperature Variation of the Top Surface Center Node 43 

2 0 Geometry of Rectangular Model 49 

2 1 Geometry of Rectangular Model 58 



ix 



LIST OF SYMBOLS 



Symbol 



Description 



a 

Al i.j.k 

A 2l.j.k 

A 3l.j.k 

B 'l.J.k 
B 2i.j.k 

Bs l.j.k 

Bi 

c 

c 1 1.j. k 

C 2l.j,k 

C 3l.j.k 

Csa 

d 

D1 i.j.k 



D 



2 i,j.k 



D 



3 i.j. k 



albedo coefficient 
coefficient associated with 

coefficient associated with 

coefficient associated with 
coefficient associated with 
coefficient associated with 

coefficient associated with 



T*j j k in Brian’s technique 
k in Brian’s technique 

kl in Brian’s technique 
T* k in Brian’s technique 
T*j* k in Brian’s technique 

T"* k in Brian’s technique 



Biot number 

specific heat of material 

coefficient associated with T*j k in Brian’s technique 
coefficient associated with T kj *, k in Brian’s technique 

coefficient associated with T"J k+1 in Brian’s technique 
solar aspect coefficient 

sum of satellite altitude and radius of earth 

coefficient associated with T* simultaneous equations in 

Brian’s technique 

coefficient associated with T** simultaneous equations in 
Brian’s technique 

coefficient associated with T n+1 simultaneous equations in 
Brian’s technique 



x 



p 

ALBEDO 

P 

COND 

P 

CONV 

^ EARTH 
P 

FLUX 

P 

GEN 

^G itjtk 

^ IN 
^ OUT 

Erad 

P 

^ SOLAR 
^ ST 

Fa 

Fe 

Fo 

h 

i 

j 

k 



rate of energy transfer into a control volume from albedo 
flux 

rate of energy transfer into a control volume from 
conduction 

rate of energy transfer into a control volume from 
convection 

rate of energy transfer into a control volume from earth 
thermal radiation 

rate of energy transfer into a control volume from imposed 
flux 

rate of energy generation 
rate of energy generation per unit volume 
rate of energy transfer into a control volume 
rate of energy transfer out of a control volume 

rate of energy transfer into a control volume from radiation 
heat transfer 

rate of energy transfer into a control volume from solar 
flux 

rate of increase of stored within a control volume 

geometric factor for albedo flux 

geometric factor for earth thermal radiation 

Fourier number 

heat transfer coefficient 

x-coordinate of node 

y-coordinate of node 

z-coordinate of node thermal conductivity of material 



xi 



p 

q" 

Ri 



R2 

T„ 

yn 




T* 



yn+ 1 
Ui 

u 2 

U 3 

u 4 

U 5 

Yi 

a 

At 

Ax 

Ay 

Az 

e 

M- 



period of revolution for earth satellite 
imposed heat flux 




ambient temperature 
temperature at previous time level n 

temperature at half-time step interval 

temperature at n + j time level in the x-direction 

temperature at n + \ time level in the y-direction 

temperature at current time level 

RlFo 

R 2 Fo 

2R x Fo 

2R2F0 

2Fo 

FoAx 2 /k 

solar absorptivity 
time increment 

distance increment in x-direction 
distance increment in y-direction 
distance increment in z-direction 
emissivity; absorptivity of earth flux 
gravitational constant 



xii 



p 

a 

<t>S 

<{>E 



mass density 

Stefan-Boltzmann constant 
solar flux 

earth thermal radiation flux 



The subscripts below have the following meaning: 



1 


left face 


2 


right face 


3 


bottom face 


4 


top face 


5 


front face 


6 


back face 


Xi 


- 2 FoaeiRiAy/k 


X 2 


-2FoG£2RiAy/k 


X 3 


-2FoG£3R2Az/k 


X4 


-2FoGE4R2Az/k 


X 5 


-2FoG£5Ax/k 


X 6 


- 2 FoG£ 6 Ax/k 


Vi 


1 

FoBijR; 


V 2 


1 

FoBi 2 R; 


V 3 


1 

FoBi 3 Rj 


v 4 


1 

FoBi . R^ 


V 5 


FoBis 


V 6 


FoBi 6 


Wi 


Ay , 

2 FoR , l (q" + 



xiii 



w 2 2FoR lir (q'; + e 2 F E (J) E +a 2 C SA 2 (t) s +a 2 aF A (J) s ) 

W 3 2FoR 2 -^-(q" + e 3 F E <J> E + oc 3 Csa 3 <1> s + a 3 aF A <J> s ) 

W 4 2FoR 2 ^(q'; + e 4 F E (J) E +a 4 CsA 4 <i) s +a 4 aF A ({) s ) 

W 5 2Fo^(q” + e 5 F e « e + cigCsA 5 4 > 5 + o 5 aF A 4> s ) 

W 6 2Fo-^(q" + £ 6 F E (| 1 E + o 6 Cs» 6 0 s + a t aF A <t>J 



XIV 



ACKNOWLEDGMENTS 



The author would like to express his deep gratitude to Assistant 
Professor Y. Joshi, his thesis advisor, for his guidance and advice in 
the endeavor. 



xv 



I. INTRODUCTION 



A. STATEMENT OF PROBLEM 

Three-dimensional transient heat conduction computations 
involving finite difference techniques can be costly and time consum- 
ing. In the explicit technique, there is a maximum time increment 
that can be used before the solution becomes unstable [Ref. 1], This 
time increment results from an upper limit on the allowable values of 
the Fourier number [Ref. 2[. The implicit technique, on the other 
hand, requires more complicated computations than the explicit 
technique due to matrix inversion operations [Ref. 2]. Even the stor- 
age capability of a modern computer can be taxed when dealing with 
storage and manipulation of very large matrices [Ref. 3], Both the 
explicit and implicit techniques can become expensive when a large 
number of nodes is required [Ref. 4]. 

For two-dimensional transients, an alternative technique is the 
alternating direction implicit (ADI) method. Calculations at each time 
step require only the handling of tridiagonal matrices [Ref. 5]. A 
direct extension of the ADI technique to three-dimensional problems, 
however, has been found to have only conditional stability [Ref. 5], A 
different algorithm by Brian is unconditionally stable and still requires 
the handling of only tridiagonal system of equations, thus avoiding the 
more complex matrix manipulation of the implicit approach [Ref. 3], 
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B. OBJECTIVES 



The objectives of this study were: 

1. To develop a FORTRAN computer program employing Brian’s 
technique for three-dimensional transient heat conduction. 

2. To validate the FORTRAN computer program using the explicit 
technique with identical geometric configuration, boundary, and 
initial conditions. Comparisons with closed-form solutions avail- 
able for a semi-infinite geometry were also made. 

3. To employ the above technique in the application of temperature 
computations involving an orbiting satellite. 

This study was also intended to be a basis for future satellite heat 
transfer analysis and computations involved with Project Orion, the 
Naval Postgraduate School’s multi-disciplinary space research 
program. 

All objectives were achieved. 
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II. DEVELOPMENT OF MODEL 



A. MODEL CONDITIONS 

The model is a three-dimensional rectangular block with geomet- 
ric dimensions expressed in Cartesian coordinates. A numerical tem- 
perature solution only permits the determination of the temperature 
distribution at discrete points or nodes [Ref. 2]. The method of finite 
differences requires that a nodal equation be written for each node in 
the block. The number of nodes in any coordinate direction is based 
on the distance increment in the same direction. 

The nodal temperature equation for each node was obtained using 
a conservation of energy equation within a control volume. Incopeia 
and DeWitt [Ref. 2] present the general form of the conservation of 
energy equation as 

^ IN _ ^OUT + ^GEN — ^ ST ( 1 ) 

The nodal temperature equations were developed using equation 1. In 
the following sections, an equation for a node on the left surface of the 
rectangular block away from corners and edges is derived for both the 
explicit and Brian techniques. Equations for all other nodes are 
derived in an identical fashion. The explicit technique equations are 
tabulated in Appendix A and the Brian technique equations are tabu- 
lated in Appendix B. 
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B. EXPLICIT TECHNIQUE 



1. Derivation of Nodal Temperature Equations 

Figure 1 shows the left surface of the rectangular block with 
associated boundary conditions. Equation 1 is rewritten assuming that 
all the energy components are into the control volume as 



where the energy generation term is assumed to be zero. 

Using finite difference notation, the energy storage term is 



The energy inflow term, E IN , is, in general, composed of convection, 
conduction, radiation, and imposed heat flux terms. For the develop- 
ment of the explicit technique, the radiation term was neglected. In 
the explicit technique, all involved nodal temperatures on the left side 
of equation 2 are evaluated at the prior time step. The convection 
term is thus 



GEN 



= E 



ST 



( 2 ) 




(3) 



E 



CONV 



= h Ax Az (T„ -T" k ) 



(4) 
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Figure 1. Left Surface of Rectangular Block Used in the Model 

The net conduction term will have contributions from all three coor- 
dinate directions and is 



_ n __ n _ _ n __n n 

e cond . k f A2 T '-^ + T ^- 2T '^ + 

COND 2 Ax Ay 



+ kAx 



Ay + T,j k+1 - 2T" j k 

2 Az 



(5) 



The flux term will be 



Ennx = q"AxAz 



( 6 ) 



Substituting equations 3 through 6 into equation 2 produces 
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*-p n , 'y ' 11 O ''p n ^P n ''P n 

k ^ Az _ ±±1 ^ ^ + kAxAz — — 

2 Ax Ay 



+ kAx 



Ay T". k . 1 +T" lJ . k , 1 - 2 T" uk 



Az 



+ kAxz(T_ -T" J k ) ( 7 ) 



Ay 

+ q"Axz = pC Ax-^-Az 



f rplH"! rp \ 

1 i.j.k ” 1 i.j.k 

At 



The temperature at the next time level, n+1, is next 
expressed in terms of temperatures at the previous time level, n. To 



accomplish this, terms in equation 7 are rearranged. 



Multiplying each term in equation 7 by 



Ax 



kAyAz 



yields 



+t ".,j.k - 2T u. k ) + (^•) 2 ( T " J .,. k -<,*) 



+ i ( f -) 2 ( T " J . k - 1 +T ,", k .,- 2T " J .0 



( 8 ) 



+ k Ax ' 0 ( T -- T " Jk ) + 5 ' 



. Ax 2 i pC 



k Ay 



kAy 2 kAt 



■^(tU -T"o. k ) 



Equation 8 is simplified using Biot and Fourier numbers, 
which are written as and 77777T. respectively. The terms 



k 

2 



pC Ax 



(#) • (f ) • and ^ 



^-^-J are identified as Ri, R2, and R 1 , respec- 
tively. Equation 8 now becomes 



6 



( 9 ) 




Rearranging terms, the temperature at the current time level is 
expressed as 



In a similar manner, the temperature equations for each node 
are derived. 

2. Restrictions on Time Step in the Explicit Technique 

In the derivation of the explicit technique equations, the 
Fourier number was expressed as 



The stability of the temperature solution is dependent on the time 
increment [Ref. 2]. Rearranging equation 10 to isolate the node of 
interest at both the current and previous time levels yields 





( 10 ) 





(ID 
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( 12 ) 



■C = Fo(T:. 1Jik + t: + 1j tk ) + 2 FoR +ltk + FoR 2 (T" j k l +T^ , k+1 ) 



Ay 



+ 2 Fob ,r;t„ + 2FoR 1 q ,/ — 

+ (l - 2Fo - 2FoR j - 2FoR 2 - 2FoBiR 



4\n 



The coefficient associated with the previous time level at the 
node of interest must be greater than or equal to zero [Ref. 2]. This 
implies 



1 -2Fol^l + Rj +R 2 + BiR; 



>0 



(13) 



Rearranging terms, it follow that, for stability. 



Fo < 



2(1 + Rj + R 2 + BiR ; 



(14) 



Limitations on Fo were derived for each different node. The maximum 
allowable Fo was chosen as the least of all these values. Since the 
Fourier number is dependent on At, Ax, and material properties, it 
follows from equation 1 1 that the maximum time increment will be 
dependent on the distance increment. As Ax decreases, the maximum 
time increment must also decrease to have equation 14 hold true 
[Ref. 1]. 
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3. Numerical Calculations 

A computer program (program EXPLICIT) was written to 
evaluate the transient temperature distribution within the block using 
the explicit technique. Figure 2 is a flow-chart description of the 
program. A complete listing of program EXPLICIT is provided in 
Appendix C. 

C. BRIAN S TECHNIQUE 

The Brian technique provides unconditional stability, allowing 
large time increments for transient three-dimensional computations 
[Ref. 3]. As mentioned earlier, the explicit and fully implicit tech- 
niques may become unsuitable for computations involving a large num- 
ber of nodes or long transient times. The explicit technique relies on 
small time increments for stability, thus increasing the amount of 
computer time required for calculations [Ref. 5], The implicit tech- 
nique requires inverting large matrices that may create problems with 
storage capability and computer time [Ref. 4]. 

Another possible technique, the alternating direction implicit 
(ADI), is applicable to both two- and three-dimensional problems. 
While the ADI method is unconditionally stable for two-dimensional 
cases, instability occurs for large time increments in three-dimen- 
sional problems [Refs. 4, 5]. 

The Brian technique modifies the ADI method to provide uncon- 
ditional stability and minimize computer time [Ref. 3]. 
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Figure 2. Flow Chart for Program EXPLICIT 
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1. Derivation of Nodal Temperature Equations 



Figure 1 shows the left surface of the rectangular block with 
associated boundary conditions. Equation 2 is again the starting point 
for nodal temperature equation derivation and is 

^ IN ^GEN “ ^ ST (2) 

where the energy generation term is again assumed to be zero. 

The nodal temperature equation derivation for Brian’s tech- 
nique uses half-time steps for each coordinate direction. This differs 
from the explicit technique, in which a whole-time increment is used 
and spatial coordinate direction is not considered. 

In Brian’s technique, the first half-time level, n + j, treats 

temperatures in the x direction as unknown with temperatures in the 
y and z directions known from the previous time level, n. Solving a 
tridiagonal system of simultaneous equations, the temperatures at the 
n + \ level are determined. 

Temperatures in the y direction are then determined using 
x-direction temperatures at the n + \ time level and z-direction tem- 
peratures at the M time level. Similarly, temperatures in the 
z-direction at the current n+1 time level are determined using n + \ 

temperatures from the x and y directions. 

The above procedure is demonstrated for a node on the left 
surface away from corners and edges. 



The energy storage term is 



Ay 

E st =P CA x— AZ 



n+— 

rp 2 rpti 

1 i,j,k “ 1 i.j.k 

At 

2 



(15) 



The energy inflow term, E IN , is composed of convection, conduction, 
radiation, and any imposed external heat fluxes. As a first step, only 
the nodal temperatures associated with heat conduction and energy 
storage in the x-direction are evaluated at the n + \ time level. The 

convection term is thus 



E C o»v=hAxAz(T.-T^ k ) 



(16) 



The net conduction term will have contributions from all three coor- 
dinate directions and is 



Ay 






-T 



n 



n+ T n +4 n + 4 

T +T -2T 

„ , . l-l.i.k T 1 WJ.k A I.j.k . . . ^ i.j+l.k * l.J.k 

E co ro = k - 2 “ Az S +kAXAZ S 



(17) 



+ kAx 



Ay + T" j k+1 ~ 2 T” j ; 



Az 



The radiation term, assuming extensive surroundings, will be 

Erad = OCAxAz(t 1 - T k j k ) (18) 

For applications involving an orbiting spacecraft, additional input flux 
terms need to be considered. These include the solar flux, the 
infrared emission from the earth, and the reflected component of the 



12 



solar flux from the earth’s atmosphere, the albedo. The solar flux term 
will be 



E souu* = Ax AzotC 



SA^S 



The reflected solar flux (albedo) term will be 



ALBEDO 



= AxAzF A a a<{) s 



The earth flux term will be 



E EARTH = AX AZF E £<|> e 



The imposed heat flux term will be 



E flux = AxAzq" 



Substituting equations (15) through (22) into equation 2 yields 



n+i n+i n+i 

#np 2 /np 2 q rnp 2 rp n rp H 

Ay 1 1 -i.j.k + 1 l+l.j.k ” Z 1 i.j.k 1 i,j+l.k “ 1 i.j.k 

k-r-Az + kAxAz 

2 Ax Ay 



+ kAx 



Ay T” j k _j +T" j k+1 -2T n 



Az 



i.j.k , . n 4 \ 

h oeAxAz^T.. - T t j k J 



+ AxAzaC <J) S AxAzF aa<t> s + AxAzF e<J) E + AxAzq" 



+ hAxAz(T„ -T"j k ) = pCAx-^-Az 



n+- 

m 2 rp n 

Ay 1 i,j.k “ 1 i.j.k 



At 

2 



(19) 



( 20 ) 



( 21 ) 



( 22 ) 



(23) 
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The temperatures at time level n + ^ are next expressed in 
terms of temperatures at the previous time level, n. To accomplish 
this, terms in equation 23 are rearranged. 

Av 

Multiplying each term in equation 23 by yields 



tIT 



* m * 

+ T 

1-1 ,j,k A 1+1 .j.k 






l.J+l.k 




1 f / ,p n p n n np ^ \ , CT0 Ax f 4 rp n 

*YAz ) V 1 u.k -1 + 1 i.j.k + i - z 1 i.j. k ; + kAy V A ~ “ A '-J k 



aC Ax 2 aF, aAx 2 F_ eAx 2 A v 2 
kAy kAy kAy kAy ^ 



(24) 



hAx 2 „ hAx 2 / T _ T n \ _ 2pCAx 2 ( n +* n 

kAy " kAy \ “ ‘j.kJ 2kAt v 1 'J- k ‘J- k 



Equation 24 is simplified using Biot and Fourier numbers, 
which are written as -j^-Ax and 0j ^f^ r2 . respectively. Equation 24 



2pCAx 



now becomes 



j( T IZ> + T"; k - 2T°**) + R, (T k - T ”, k ) 

+ TR 2 (T n LJ . k _, +T", k+1 - 2T”, k ) + f-R,Ay(Ti - T k ; k ) 



aC ^ R i A y aF A aR,Ay F E eR,Ay 



<>s + 



4>s + 



■Of 



(25) 



R 



1 A y „ 2 / n \ 1 ( n+ 2 n 

k ^ + B i R i (T„ -T, 0 k ) = 2 Fo l T ‘+ k “ T i.j.kJ 
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Rearranging terms, the temperatures at the first half-time step level 
in the x direction are 



n+ r 



- F ° T ,-, 0 .„ + (1 + 2Fo)T m ’ - FoT 1Jk = T- v 



j+l.k 



+ 2FoR,(t " j+1 k "T", k ) + FoR 2 (T" kk ., + 4, - 2T kJ k ) 



aC 



aF. a 

A 



F e 

1 E C 



+ 2 Fo (}) s + 2Fo R, Ay-^- 4) s + 2 Fo R, Ay-^ (J) E 



(26) 



+ 2Fo RjAy + ^ + 2FoBi R* (T„ -T” k ) 



The next step is to express the temperatures in the y-direc- 
tion at time level n + y using x-direction temperatures determined at 
time level n + 7 from equation 26 and z-direction temperatures at 
time level n. For greater clarity, the x-direction temperatures at time 
level n + j, T 2 , are rewritten as T*. Equation 26 now becomes 



-FoT*, Jk + (1 + 2Fo)T,“ k - FoT* +1 j k =T” J+1 k 

+ 2FoR,(T" j+1 k -Tj'j k ) + FoR 2 (T kJ + T kJ k+1 - 2T kJ k ) 

ne/ , n 4 \ ctCL. aF. a (071 

+ 2Fo R,Ay^(T. - T tJ J + 2FO R.Ay-jf- $ s + 2FO R,4y— j^-$ s 



F E e 



+ 2 Fo R , Ay 0 E + 2 Fo R , Ay -^ + 2 FoBi R ’, ( T . - T" , k ) 



The derivation for the y-direction nodal temperature equation is simi- 
lar to the x-direction derivation. The y-direction equation for the same 
node is 
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n+ 7 



n+ T 
2 



(1 + 2FoR,)T,^ -2F0R.T.J,, , k = T, J k 

+ F 0 (T* 1J . k +T *,j k “ 2T t * k ) +FoR 2 (T kJk _j +T kJk+1 -2T kJk ) 
+ 2Fo R.Ay^T 4 . - t;‘ k ) + 2Fo R.Ay^ t> s 



cxF a F„ e 

+ 2 Fo Rj Ay— <J) S + 2Fo RjAy — (]) E 



(28) 



+ 2FoR 1 Ay^ + 2FoBi R^T.-T^) 

1 ® 

The y-direction temperatures at time level n + ^ . T , are rewritten 
as T**. Equation 28 now becomes 

(1 + 2FoR ,)T” -2FoR l T”. k = T” k(+ Fo(T* u „ + T;„ lk -2T* JJt ) 

+ FoR 2 (T, n Jk . l +T" Jk „ -2T, n Jk ) +2Fo R,Ay ^( t! - T^ k ) 



ccC 



aF A a 



F e e 



va.v^q* ua a <-*. *■ p tx 

+ 2Fo Rj Ay — ^ — <|> s + 2Fo R, Ay — ^ — <{) s + 2 Fo R ,Ay — <I> E 



(29) 



+ 2Fo R, Ay y + 2FoBi R, (T. -T ” ik ) 



The expression for the z-direction temperatures at time level 
n +1 is dependent on the x-direction and y-direction temperatures at 
time level n + 7 , T* and T**, respectively. The z-direction nodal tem- 
perature equation is 
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- for 2 t-L + (i + 2 For 2 )t;;;; - for 2 t;;;; = 



n+1 



,n+l 



+Fo ( T ,t 10 . k + T* 1J .v -2T*. k )+2FoR 1 (T*; i „ - T” ) 



aC. 



+ 2FoR,Ay^(Tl - T" ) + 2Fo R,Ay-jfH s 



aF A a 



F E e 



+ 2Fo RjAy k 4> s + 2FoR 1 Ay^-({) E 



( 30 ) 



+ 2Fo R,Ay-^- + 2FoBiR J (T. - T* * ) 



In a similar manner, the temperature equations for each node 
are derived. We note that during traverse along each coordinate direc- 
tion in equations 26, 28, and 30, a tridiagonal matrix results for the 
solution of temperature. 

2. Solution of the Nodal Temperature Equations 

The x-direction nodal temperature equations for the rectan- 
gular block are of the form 

Aq.j.jr-,,, + ® h.j.k'F *j.k + ^' 1 i.j.k^'i+l.j.k = ^ 1 l.j. k (31) 

where the coefficients Ai t jk , .... D i 1Jk are known from the previous 
time step. 

Each row of T* equations parallel to the x-axis is a tridiagonal 
matrix of coefficients Ai ljk , Bi ljk , and Ci ljk . It was solved using the 

standard Tridiagonal Matrix Algorithm (TDMA) [Ref. 5], After each row 
of T* equations is solved, the y-direction equations will be solved [Ref. 
3], These are of the type 
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(32) 



A 2 ,, k Tr j !, u 



+ B 



9 T* * 

2 i.j.k i.j.k 



+ c 



9 T* * 

2 t.J.k i.J+l.k 



= d 2 



i.j.k 



As in the x-direction case, each row of T** equations parallel to the 
y-axis is solved using a tridiagonal matrix solver. The z-direction equa- 
tions are of the form 



A3 i.j.k T ”.J.k-l + B 3 i.j.k T u.k + C3 i,j.k T ij.k + l - D3 i.j,k 



,n+l 



n + 1 



(33) 



The solution of the T n+1 equations is accomplished in a similar man- 
ner as used for both the x- and y-directions [Ref. 3]. 

3. Numerical Calculations 

A computer program (program (BRIAN) was written to evalu- 
ate the transient temperature distribution within the block using the 
Brian technique outlined in the previous sections. Figure 3 is a flow 
chart description of the program. A complete listing of program 
BRIAN is provided in Appendix D. 
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Figure 3. Flow Chart for Program BRIAN 



19 



III. VALIDATION 



A. SELECTED GEOMETRY AND CONDITIONS 

Both the Brian’s technique and explicit technique programs 
required validation to ensure accuracy in the temperature solutions. 
Identical geometric and boundary conditions were used for these 
techniques in the validations. The left face had a constant surface heat 
flux with all other surfaces treated as adiabatic. Figure 4 shows the 
geometric and boundary conditions for the validation model. Table 1 
lists the parameters used in the validation runs. 

The Brian’s technique and explicit technique programs were 
written to allow the user to conduct a validation test. The test involved 
taking the temperature differences at identical nodes in each program 
and then comparing each node’s temperature difference on a graph. 

The finite difference solutions for the rectangular block were also 
compared with an analytical solution for a semi-infinite solid exposed 
to a uniform flux at its surface. The closed-form solution would be a 
reasonable approximation to the actual solution for small times when 
interior nodes are not influenced by the lateral boundary surfaces of 
the rectangular block [Ref. 2j. 

Incropera and DeWitt [Ref. 2] tabulated the closed-form solution 
for the case of constant surface heat flux in terms of complementary 
error functions. These were evaluated in the present study based on 
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Figure 4. Rectangular Model Showing Boundary Conditions 

Used in Validation Runs 

TABLE 1 

PARAMETERS USED IN VALIDATION RUNS 



Left face flux 


1,500 W/m 2 


Other faces 


0 W/m 2 (adiabatic) 


Material density 


2,770 kg/m 3 


Thermal conductivity 


177 W/m-K 


Specific Heat 


875 J/kg-K 


Initial Temperature 


275K 


Ambient Temperature 


295 K 


Distances from face 


1.0, 0.4, 0.0 m 



the polynomial approximations provided by Abramowitz and Stegun 
[Ref. 5]. 
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B. NUMERICAL COMPUTATIONS 

A computer program was written to evaluate the closed-form 
solution for a semi-infinite solid with a constant surface heat flux. The 
fifth-order polynomial approximations to the complementary error 
function have an error less than 1.5 x 10' 7 [Ref. 5]. 

Figure 5 is a flow-chart description of the closed-form solution 
program VALID. The parameters entered are identical to those used in 
programs EXPLICIT and BRIAN. A complete listing of program VALID 
is provided in Appendix E. Table 1 lists the parameters used for the 
validation runs. 

The validation period for programs VALID, EXPLICIT, and BRIAN 
was 3,600 seconds (one hour). The distance inside the semi-infinite 
solid was based on pre-selected nodes in the rectangular model. Three 
distances were used in the semi-infinite solid, corresponding to three 
selected nodes in the rectangular block. 

In the first validation run, each coordinate direction had 11 
nodes. This corresponded to an explicit technique time increment of 
17.117 seconds, which was also selected as the time increment for 
programs BRIAN and VALID. In the second validation run, the number 
of nodes was increased to 21 in each coordinate direction. The time 
increment used in this run was 4.279 seconds. A third validation run 
was conducted using 1 1 nodes in each coordinate direction. The 
explicit technique and program VALID used a time increment of 
17.117 seconds. The time increment for program BRIAN was 



22 




Figure 5. Flow Chart for Program VALID 
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increased to 120 seconds. The increased time increment was used to 
determine whether the program BRIAN temperature differences 
would differ significantly from those previously obtained with the 
17.117-second time increment. 

C. VALIDATION COMPUTATIONS 

Graphical representations of the temperature differences experi- 
enced over the validation period are presented in Figures 6 through 
14. 

At a distance of one meter into the block (Figures 6 and 7), the 
temperature differences calculated in EXPLICIT and BRIAN diverge 
significantly from VALID after a period of 616 seconds. This is to be 
expected as the effects of the constant surface heat flux will provide 
higher temperatures in the rectangular block than in the semi-infinite 
solid. The rectangular block can thus be modeled as a semi-infinite 
solid for only about 10 minutes at a distance of one meter. 

At a distance of 0.4 meters into the solid and on the surface 
(Figures 8 through 11), both EXPLICIT AND BRIAN temperature dif- 
ferences were in close agreement with the closed-form temperature 
differences for almost the entire validation period of one hour. 

Increasing the number of nodes and, consequently, decreasing 
the time increment had virtually no effect on temperature differences, 
as can be seen in Figures 6 through 11. 

Increasing the program BRIAN time increment to 120 seconds 
did not degrade the accuracy of the solution. In the three cases 
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Figure 6. Temperature Difference Versus Time for Time Increment 
of 17.117 Seconds with Distance of 1 Meter Into Block 
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Figure 7. Temperature Difference Versus Time for Time Increment 
of 4.279 Seconds with Distance of 1 Meter Into Block 
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Figure 8. Temperature Difference Versus Time for Time Increment 
of 17.117 Seconds with Distance of 0.4 Meters Into Block 
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Figure 9. Temperature Difference Versus Time for Time Increment 
of 4.279 Seconds with Distance of 0.4 Meters Into Block 
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Figure 10. Temperature Difference Versus Time for Time 
Increment of 17.117 Seconds on Left Surface of Block 
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Figure 11. Temperature Difference Versus Time for Time 
Increment of 4.279 Seconds on Left Surface of Block 
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Figure 12. Temperature Difference Versus Time for Time Increments 
of 17.117 and 120 Seconds with Distance of 1 Meter Into Block 
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Figure 13. Temperature Difference Versus Time for Time Increments 
of 17.117 and 120 Seconds with Distance of 0.4 Meters Into Block 
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Figure 14. Temperature Difference Versus Time for Time 
Increments of 17.117 and 120 Seconds on Left Surface of Block 
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involving distances of 1.0 and 0.4 meters into the surface and on the 
surface, the temperature differences did not vary by more than 0.045 
percent from those obtained using a 17.117-second time increment. 
Figures 12 through 14 show the 120-second time increment plot of 
temperature difference versus time being in close agreement to the 
corresponding 17.117-second time increment plots. 

Using a time increment of 17.117 seconds, the maximum error 
was found to be within 0.019 percent for an 11 x 11 x 11 grid. The 
maximum error was calculated using 

% ERROR = Tcompu ™ Iexact. x 100% 

* EXACT (34) 

where Tcomputed is the Brian’s technique computed temperature and 
Texact is the closed-form temperature. Sampling points were taken 
at 600-second intervals. 

In a similar manner, the maximum error was calculated for a time 
increment of 4.279 seconds. The maximum error was found to be 
0.021 percent for a 21 x 21 x 21 grid. 

A time increment of 120 seconds was above the maximum allow- 
able time step from stability considerations for the explicit technique. 
Brian’s technique thus allowed time steps above the limitation of the 
explicit technique with very little compromise in solution accuracy. 

Based on the validation test results, program BRIAN can be con- 
sidered correct and validated for further use. 
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IV. APPLICATIONS 



A. GENERAL 

The Brian technique was used to solve for the temperature distri- 
bution in an application example. The example involves the transient 
thermal response of a component on board an orbiting satellite. 

B. ORBITING SATELLITE COMPONENT 
1. Geometry and Initial Conditions 



The transient thermal response of an orbiting spacecraft 



component will be investigated. The spacecraft component is in a non- 
geosynchronous circular orbit at an altitude of 1,609.364 kilometers 
(1,000 miles). The equation for the orbital period is listed in Agrawal 
[Ref. 7] as 



where p is the gravitational constant of the earth and d is the sum of 
the satellite altitude and radius of the earth. Eisele and Nichols [Ref. 8] 
list the values of p and earth radius as 398,603.2 km 2 /s 2 and 
6,378.165 kilometers, respectively. 

The period of revolution in this case is 7,104.4 seconds or 
118.40 minutes. For a satellite in a circular orbit at an altitude of 100 
miles, Stevenson and Grafton [Ref. 9] listed the sun exposure time as 
72 percent of the orbital time. During one revolution around the earth, 




(35) 
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the satellite will have 85.25 minutes of sun exposure and spend 33.15 
minutes in the earth’s shadow. For the example application, the satel- 
lite is assumed to start the orbital period at the instant it becomes 
exposed to the sun. 

Single spin stabilization is assumed, thus allowing the 
spacecraft component to spin about a principal moment of inertia axis 
[Ref. 7]. The spacecraft component spins about the z-axis, as depicted 
in Figure 15. With single spin stabilization, it is assumed that the 
component analyzed here will have all four sides continuously exposed 
to the sun with the top and bottom surfaces not receiving solar flux. 
The fraction of solar flux that is incident on the satellite surface is 
usually defined as the solar aspect coefficient [Ref. 7). Thus, for this 
example, the solar aspect coefficients are assumed to be 1.0 for the 
side surfaces and 0.0 for the top and bottom surfaces. 

Stevenson and Grafton [Ref. 9] compiled the geometric fac- 
tors for both earth radiation and reflected solar radiation (albedo) on a 
flat plate. The geometric factor for earth radiation depends on the sat- 
ellite’s altitude and the altitude angle between the normal to the rect- 
angular component and the vertical to the satellite from the earth [Ref. 
9]. Figure 16 depicts the geometry involved. The geometric factor for a 
rectangular plate at an altitude of 1,000 miles (1,609.364 kilometers) 
and altitude angle of 60 degrees is listed as 0.3174 [Ref. 9]. 
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Figure 15. Component on a Satellite in Circular Orbit 
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ALTITUDE 




Figure 16. Geometry Used in Determining Geometric Factor 
for Earth Thermal Radiation 

The geometric factor for reflected solar radiation is depen- 
dent on the satellite’s altitude, the altitude angle, the sun angle, and 
the angle of axis rotation [Ref. 7], Figure 17 depicts the geometry 
involved. The geometric factor is 0.0456 for a satellite at 1,000 miles 
altitude, altitude angle of 60 degrees, sun angle of 90 degrees, and 
rotation angle of 0 degrees [Ref. 9]. 
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Figure 17. Geometry Used in Determining Geometric Factor 
for Solar Reflected Radiation 
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TABLE 2 



PARAMETERS USED IN ORBITING SATELLITE APPLICATION 



Dimensions 

Material 

• density 

• thermal conductivity 

• specific heat 

Exterior Coating 

• solar absorptivity 

• emissivity 

Imposed Flux 

• bottom surface 

Orbital Data 

• altitude 

• period of rotation 

• exposure to sun 

• exposure to darkness 

Solar Aspect Coefficients 

• side surfaces 

• top and bottom surfaces 

External Flux 

• solar flux 

• earth radiation 

• albedo coefficient 

Geometric Factors 

• earth radiation 

• albedo flux 

Initial Data 

• heat transfer coefficients 

• ambient temperature 

• initial nodal temperature 



0.15 m x 0.10 m x 0.015 m 

aluminum alloy 2024-T6 
2,770 kg/m 3 
186 W/m-K 
1,000 J/kg-K 

black paint 
0.97 
0.92 



50,000 W/m 2 

1,609.364 km (1,000 mi) 
118.40 min 
85.25 min 
33.15 min 

1.0 

0.0 



1,353 W/m 2 
237 W/m 2 
0.3 



0.3174 

0.0456 

0.0 W/m 2 -K 
0 K 
500 K 
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Agrawal (Ref. 7] listed the mean annual solar flux, mean 
annual earth radiation, and mean annual value of albedo coefficient as 
1,353 W/m 2 , 237 W/m 2 , and 0.30, respectively. These value are used 
in the example. 

The material coating on the component's exterior determines 
the solar absorptivity and emissivity. The component is assumed to be 
coated with black paint, thus providing a solar absorptivity of 0.97 and 
an emissivity of 0.92 [Ref. 10]. 

The component is rectangular with dimensions of 0.25 m x 
0.10 m x 0.015 m and is made of aluminum alloy 2024-T6. Incropeia 
and DeWitt (Ref. 2] list the values for density, thermal conductivity, 
and specific heat. 

The component is mounted on the exterior of the satellite. 
The bottom surface is attached to the satellite with all other surfaces 
exposed to the space environment. The flux on the bottom surface is 
assumed to be 50,000 W/m 2 . 

Table 2 lists the properties and other parameters used in the 

example. 

2. Numerical Calculations 

A time increment of 120 seconds (2 minutes) was used in 
program BRIAN. The program covered a period of 120,000 seconds 
(33.33 hours), during which the spacecraft component conducted 
16.9 orbits. The temperatures of top surface node were plotted versus 
time, as shown in Figures 18 and 19. The temperature versus time 
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Figure 18. Temperature Variation of the Top Surface Center Node 
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Figure 19. Temperature Variation of the Top Surface Center Node 
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plots in both figures demonstrate an imperfect sinusoid. Steady peri- 
odic response occurs following orbit. The peak temperatures during 
the third orbit vary from the second orbit peak temperatures by 0.27 
percent. The largest temperature difference in successive peak tem- 
peratures is 0.27 percent and occurs between the second and third 
orbits. 

The oscillations in temperature distribution may be attributed 
to overlaps in the sunlight and eclipse periods. The sunlight exposure 
time for the orbit is 85.25 minutes. A time increment of 120 seconds 
(2 minutes) would allow solar flux to be seen by the satellite compo- 
nent for an additional 1.25 minutes on the first orbit. Thus, 1.25 min- 
utes of the first eclipse period would be missed at the beginning of the 
eclipse period. At the end of the eclipse period, the component is in 
darkness an extra 1.60 minutes. These overlaps occur throughout the 
running of the program with the maximum overlap being 2 minutes. 
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V. RESULTS 



Brian’s technique was validated for four of the six rectangular 
surfaces. For cases involving flux or convective boundary conditions on 
the front and back surfaces of the model, temperature distributions 
did not follow the pattern obtained in previous validation cases. This 
anomaly requires further investigation. 
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VI. CONCLUSIONS 



Brian’s technique was validated using the explicit technique and 
the closed form solution for a semi-infinite solid. For the same time 
increment, the explicit technique uses six times less CPU time. How- 
ever, as the time increment was increased, the CPU time used in 
Brian’s technique became less than that used in the explicit technique 
for time steps larger than 110 seconds. Brian’s technique used time 
increments that would lead to instability if used in the explicit 
technique. 



46 



VII. RECOMMENDATIONS 



While developing the computer code for Brian’s technique, mate- 
rial properties were assumed to be independent of temperature. A 
more precise code would make thermal conductivity a function of 
temperature. 

Another improvement for the computer code when used for 
spacecraft thermal analysis computations involves the application of 
orbital mechanics. This would permit the use of non-circular orbits in 
the code. 

This study was restricted to the use of Cartesian coordinates in 
Brian’s technique. Further study is required in applying Brian’s tech- 
nique in cylindrical and spherical coordinates. 
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APPENDIX A 



EXPLICIT TECHNIQUE NODE EQUATIONS 



A. INTERIOR NODE EQUATION 



T m „ = Tt, k + Fo(t; i j k + T” , 4 k - 2T" J k ) 
+ R,Fo(T: J ., k+ T: Jt ,. k -2T: Jk ) 

+ R 2 Fo(T" jk _, +T kjktl - 2T" Jk ) 



B. SURFACE NODE EQUATIONS 
1. Left Face 



■c, = + fo(t” u k + t“ Uk - 2x: J . k ) 

+ 2F0R 1 (T" J „ k -T" Jk ) 

+ FoR 2 (T" Jk _, +T” Jk+) ~2T kjk ) 

1 

+ 2Foq"R j(T_ -T"j k ) 



2. Right Face 



T "i = T M. k + Fo(T“ ,, k + T”„ 4 - 2 T"j k ) 
+ 2FoR,(T° J _, k ~T"j k ) 

+ F0R 2 (T; ik _ 1+ T". k .,-2T: j . k ) 
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Figure 29. Geometry of Rectangular Model 
Bottom Face 

•C = t m.„ + + - 2T " )k ) 

+ F°Ri( T "j_i. k + T" Jtl k “ 2T" J k ) 

+ 2FoR 2 (T" J ktl -T” J k ) 



Top Face 



T 



n + 1 
l.j.k 



= T;. ) , k+ Fo(T:_ 1J . k + T;, J , k -2T: )t ) 
+ FoRjCt^j k +T kJkl k - 2T kJ k ) 

+ 2FoR 2 (T" Jk . 1 -T" Jk ) 



RIGHT 
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5. Front Face 



Ti;L=T; k+ 2Fo( T ;. 1 , k - T ; J ,) 

+ F0R 1 (T" J .,, k +T” J „„-2T: jll ) 

+for 2 (t; ) . i , -1 ^i.J.k+l 2T j j k) 

6. Back Face 

T M> T u k + 2 F K T ^^- T ”j.O 

+ FoR j(T ( j_j k + T, J+1 k - 2T (J k ) 

+ f 0 r 2 (t: j ,. 1+ t: j , +1 

2T t.j.k) 

C. PERIMETER NODE EQUATIONS 
1. Perimeter 1 

T;; 1 k = T;., k+ 2Fo(T; +1J . k -T; jk ) 

+ FoRj^T^j k +t” j+1 k - 2T” J k ) 

+ 2 FoR 2 (T "j k _j - T k ) 
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2. Perimeter 2 



T "j*k = T l!j.k + 2Fo(T"_, Jk -T°j k ) 

+ F0R,( T "j-..k +T " J .,.k- 2T ”,.k) 
+ 2FoR 2 (T” Jk „,-T" Jk ) 

3. Perimeter 3 

T:.;;=T;, k+ 2Fo(x:., J ,-T:, k ) 

+ FoR,(T 1J _, k +T 1J<lk -2T tJk ) 
+ 2FoR 2 (T" Jktl “T" Jk ) 

4. Perimeter 4 

T M*I = T U.k + 2Fo(T:., Jk -T; >k ) 

+ for 1 (t: j . 1 , + t: ji1 ,-2t; 1 ) 

+ 2FoR 2 (T” Jk4l -T” Jk ) 

5. Perimeter 5 

T ”;i,= T y.k+2Fo(T:. 1J . k -T: Jk ) 

+ 2FoR,(t" j _, k “ T kJ k ) 

+ F ° R J ( T y.k-, +T y.k„- 2T M.k) 
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6. Perimeter 6 



T Mk = T^ k+ 2F 0( T ; iJI ,-T" J . k ) 

+ 2 FoR , ( T k — T k ) 



+ FoR 



( T n +T n — 2T n ^ 

2^ 1 l.J.k-1 T 1 i.J.k+1 ^ 1 tj.k ) 



7. Perimeter 7 



T Uk = T " J .k + 2F°(Tw. 1 . l -T:. J . k ) 

+ 2FoR 1 (T”., k -T” Jk ) 

+ FoR 2 (T; jk . 1+ T" kkl -2T; jk ) 

1_ 

+ 2FoBiR * (T. - T”j k ) + 2Foq"R, 



Ay 



8. Perimeter 8 



^"j'k = T t!j.k + 2Fo(T" +1 j k “T", k ) 

+ 2F0 R ,(T”, lt -T" Jk ) 

+ F ° R 2( T M.k-, + ^i.J.k+1 ^^'l.J.k) 

+ 2FoBiR*(T.-T" J>k ) + 2Foq"R 1 J Y- 
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9. Perimeter 9 



T ”j*.k - T^k + + T"„ uk - 2T; j k ) 

+ 2FoR,(T" J . 1 . k -T" j . k ) 

+ 2FoR 2 (T" Jk+1 -T" Jk ) 



10. Perimeter 10 



t m . = T ;,.k + Focn.j.k + 'Co.k - 2T : 0 .k) 
+ 2FoRj(T, k — T 1Jk ) 

+ 2FoR a (T" J k _, -T" Jk ) 



11. Perimeter 11 



t:;! = t:,, + f 0 (t:. 1j k + - 2 t; j .„) 

+ 2FoR,(T” J . lk -T” Jk ) 

+ 2FoR 2 (T” Jk .,-T” Jk ) 

+ 2FoBiR *(T_ - T” k ) + 2Foq"R,-^- 
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12. Perimeter 12 



T"' k = T" J k + Fo(t;. 1 0 k + T" tl j k - 2T” k ) 

+ 2FoR,(T" Jllk -T" Jk ) 

+ 2FoR,(T" J k „ -T" J k ) 

+ 2F0B1R *(T_- T" k| k ) + 2Foq"R,^- 

D. CORNER NODE EQUATIONS 

1. Comer 1 

T” , k = T". k+ 2Fo(T; 41J . k -T: j . l ) 

+ 2FoR 1 (T" J , lk -T” Jt ) 

+ 2F0R J (T 1 " Jk . 1 -T" Jk ) 

Ay 

+ 2Foq"R 1 -jf 

2. Corner 2 

T” , k = T: jk+ 2Fo(T:. 1J . k -T" Jk ) 

+ 2FoR 1 (T" tlk -T" l ) 

+ 2FoR 2 (T," Jk . 1 -T” Jk ) 

Ay 

+ 2Foq'R | - j f 
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3. Comer 3 



T"L=T" tJ , + 2Fo( T ;.,, k -T:. J ,) 
+ 2FoR 1 (T I Jtl k - T 1J k J 
+ 2FoR 2 (T" Jk _ 1 -T" Jk ) 

Ay 

+ 2Foq'R 1 - j i 



4. Comer 4 



T “j i, = T" J k + 2F-0(T7_, j k - T” k k ) 

+ 2F0R 1 (T^ lk -T 1 " Jk ) 

+ 2F0R 2 (T” k .,-T” k ) 

Ay 

+ 2FoqR,-^- 



5. Comer 5 



T;;L=T:, k+ 2Fo(T;.,, k -T:, k ) 
+ 2F0R,(T" J . lk -T" Jk ) 

+ 2FoR 2 (T" J k+J - T" J k ) 



+ 2FoBiR ^(T„- T kJ k ) 
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6. Comer 6 



T;;L = T" l) . k + 2Fo(T;., J . k -T;, k ) 

+ 2FoR 1 (T” H . k -T”. k ) 

+ 2 FoR 2 ( T "j k+I - T ", k ) 

1 

+ 2FoBiR^(T.-T; jk ) 

7. Comer 7 

T n v ;L= T ” y . 1 . +2 F' o ( T ?«* k - T " ik ) 

+ 2FoR,(t" j _, k -T" J k ) 

+ 2FoR 2 (T; jk . 1 -T" Jk ) 

+ 2FoBiR j(T_- T kJ k ) 

8. Comer 8 

T ;;L=T; Jk+ 2Fo(T:. 1Jk -T; Jk ) 
+ 2FoR,(T” Hk -T" Jk ) 

+ 2FoR 2 (T" Jk .,-T" Jk ) 

+ 2FoBiR’(T.-T; jk ) 



56 



APPENDIX B 



BRIAN’S TECHNIQUE NODE EQUATIONS 

A. X-DIRECTION EQUATIONS 

1. Comer 1 

(l + 2Fo)T* k -2FoT* Jk =T" jk + U 3 (t" j+] k - T kJ k ) 
+ U 1 (T" Jk „-T" J . k )+V 1 (T. 1 -T^ k ) + V 3 (T.3-T; jk 
+ V 6 (T_ 6 -T" J k ) +W, +W 3 + W 6 + X,(T"’ k -T 4 .,) 

+ X,(Tt, k -Ti 3 ) + X 6 (t:;„ -Tt 6 ) + Y,E 0|Jl 

2. Perimeter 12 

"FoTf-u.k + (l+2Fo)T*, k -FoT„ 1Jk = 

T: J . k + u 3 (T; J4 , k -T; J . k ) +u 4 (T; Jktl -T" Jk ) 

+ V](T-i “ T"j k ) +V 3 (T. 3 -T," Jk ) +W, + w 3 
+ X ,( T M* k -T-,) +X 3( T i -O + Y,E ClJ k 
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LEFT 



3. 




RIGHT 



Figure 30. Geometry of Rectangular Model 
Corner 2 

-2F0T;_ 1Jk + (l+2Fo)T* t = T kJ k + U 3 (t" j+1 „ -T" Jk ) 
+ U 4 (T" Jk . 1 -T" Jk ) +V,(T., -T" kJk ) +V 3 (T. 3 -T” )k ) 
+ F 5 ( T . 5 - T ; j k)+W, + W 3 + W 5 + X 1 (Tr ) 4 k -Tl 1 ) 

+x 3« t - T y +x s (T;j. -t 4 . 5 ) +y,E 0 ,. ) . k 
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4. Perimeter 8 



(l + 2Fo)T* k - 2FoT* 4I j k =T",, k + U 2 (T" Jk . 1 + T” JW1 - 2T" Jk 
+ -Ty +V,(T., -T" k ) +V 6 (T. 6 -T". k ) 

+ W, +W 6 + X,(T; j i k -Tt.) +X 6 (T:j‘ k -Tl 6 ) + Y,E GlJk 

5. Left Face 

-FoTr. 1J k +(l+2Fo)T* k - FoT* +1 kk = T"j k 
+ U 2 (T” , k _. +T" kll - 2T kj k ) + U 3 (T" kl k -T" J k ) 

+ V,(T., - T kJ k ) + W, +X I (T;‘ k -Tt.) +Y,E 0kJk 

6. Perimeter 7 

-2 FoT j k + (1 + 2 Fo)T kJ k =T; j k 

+ U a (T^ k ., +T^ kkl - 2T; j k ) + U 3 (T" j+1 k -T”j k ) 

+ V .( T -. - T "j.O +V 5 (T.5- T " lk ) +W,+W S 

+ x ,«‘ k -T-,) + X 5 ( T " J '! k - T -s) +Y i E o,j. k 
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7. Corner 3 



(1 + 2Fo)T*j k - 2FoT* 41 J k = T" j k + U 3 (T" J+1 k - T kj k ) 
+ U 1 (T," J . lt . | -T" Jk ) + V 1 (T. 1 -T: jk ) 

+ V 4 (T. 4 -T n tJ . k )+V,(T. 8 -T" , k ) + W j ■+• 



+ x.f t” 4 . 


-tI.) +x,(t 


1\ i.j.k 


1/ 4 V 


/ n ^ 




+ X 6 (T,j. k 


- T~ 6 ) + YjE g 



8 . 



9. 



Perimeter 11 

-FoT*_j jk + (1 +2Fo)T* k -FoT* +l j k = T" J>k 
+ U 3 ( T V k -T; kk ) +U 4 (T; jk ., -T” Jk ) +V,(T., -T"j k ) 
+ V 4 (T- 4 -T” „) + W, +W k + X,(T” J 4 k -Ti,) 
+x 4 (Tr. J 4 k -T; 4 )+Y 1 E GlJk 

Comer 4 

-2F0T;. 1J k + (l+2Fo)T* k = T kj k + U 3 (t" j+1 k - T kJ k ) 
+ U 4(T",. k _, -T"j k ) +V,(T., -T l >k ) +V 4 (T. 4 -T” Jk ) 

+ V 5( T -5- T M k ) +W ,+ W 4 +W 5 +X ,( T "j 4 k - Ti .) 

+ X 4(T" J 4 k -T: 4 )+X s (T: j 4 k -Tl 5 ) + Y,E 0 , Jk 
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10. Perimeter 3 



(l + 2 Fo)T* k - 2 FoT*, Jk =T" J . k + U 1 (T". | . k +T" lJ „. k -2T^ t 
+ U .( T ",. k ., - T ",. k ) + V 3 ( T -3 -T”. k ) +V 6 (T. 6 -T” Jk ) 

+ W 3 + W 6 + X 3 (T; j 4 k -TI 3 ) +X 6 (T" 4 k -t! 6 ) + Y,Eo, Jk 

11. Bottom Face 

- FoT r-,.j. k +(l+2Fo)T*. k - FoT* Iik = T^„ + 

U .( T ".,-..u + Ty.,. k - 2T” y . k ) + U 4 (T'[ j , k „ - T kJ k ) 

+ V 3 (T - 3 -T" . k ) + W 3 +X 3 (T^ k -Ti 3 ) + Y,E GkJ k 

12. Perimeter 4 

-2 FoT + (1 + 2Fo)T* Jk =T; jk + 

U ,(Ty_,. k + T” kl . k - 2 T ;, k ) + U 4 (T kJ ktl - T" , k ) 

+ V 3 (T . 3 -T kJk ) +V 5 (T. 5 -T” Jk ) +W 3 + W 5 

+ x s( T ^ - T -a) +X 5 ( T i - T - 5 ) +Y,Eo,., k 
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13. Back Face 



<l + 2Fo)T* k -2FoT* 1Jk =T", k+ U 1 (T”. lk + T”„ k -2T" k ) 
+ U 2 (T". k _, + T:, k .,-2T”,) +V 6 (T. 6 -T" Jk ) + W 6 
+ X 6 (T: j 4 k -T 4 . 6 ) + Y 1 E 0 ,. J . k 

14. Interior 

-FoT*_, j k +(l+2Fo)T* k - FoT* +1 j k 

+u ,( T :,-.. k +T u., k - 2T «.o 

+ U 2 ( T 1 .j.k-1 +T IJ.k.l “ 2T l.j.k) + Y 1 E 0 |.j. k 

15. Front Face 

-2 FoT r- J.k + « + 2 Fo)T * k = T” kJ k + U,(T kj l k + T” j+ , k - 2Tl k ) 
+ U 2 ( T : j . k - 1 +T :.,.k., -2T; jk ) +V 5 (T. s -T" Jk ) 

+ w 5 +X 5 (Tr j 1 k -T; 5 ) + Y 1 E GlJk 

16. Perimeter 1 

(1 + 2Fo)T f * k - 2FoT* tl Jk =T” J . k + U,(T“. 1 , k + T”. 1 . k -2T" y . k ) 
+ U 4 (T°. k ., -T" k ) +V 4 (T. 4 -T" k ) +V 6 (T. 6 -T” Jk ) 

+ W 4 +W 6 +X 4 (T" 4 k -Ti 4 ) +X 6 (T” 4 k -Tl 6 ) +Y,E 0lJk 
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17. Top Face 



-Fon, Jk + (1 + 2Fo)T * k -FoTf„ )k =T" J k 
+ U,(T: j _,. k + T" 41 . k -2T; )k ) +U 4 (T” k _, -T; jk ) 

+ V 4 (T. 4 -T“. k ) +W 4 + X 4 (T"" k -Ti 4 ) + Y.Eo,.,., 



18. Perimeter 2 



-2FoT;., ; k + (1 + 2Fo)T * k = T kJ k + U^T^, k +T^, „ - 2T^ k ) 
+ U 4 (T” )k _, -T” Jk ) +V 4 (T.„ - T kj k ) +V 5 (T. 5 -T" )k ) 

+ W 4 + W 5 + X 4 (T; 4 k -Tt 4 ) +X 5 (T” ) 4 k -Tl 5 ) +Y,E 0 ,. lk 



19. Comer 5 



(l + 2Fo)T* k -2FoTf, ]Jk =T" Jk+ U 3 (T; j . lk -T; jk ) 
+ U 4 (T," ) . k41 -T” k ) + V 2 (T- ! -T". k ) + V 3 (T-3 -T" .„) 
+ ve(T. 6 -T; Jk ) + w 2+ w 3 + w 6+ x 2 (T; 4 k -Ti 2 ) 

+ X 3«„ - t4 - 3 ) +X 6(T:‘ k -Tl 6 ) +Y,E 0lJ k 
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20. Perimeter 9 



-FoT*_ 1)k + (1 + 2Fo)T * k -F 0 T*. 1Jk =T° J , k + l^T"^ - T” tJ .„) 
+ U.(Tr j . fcrt -T", k ) +V 2 (T^ - T kJ k ) +V 3 (T. 3 -T” k ) 

+ W 2 + W 3 +X 2 (t ; j 4 |1 -Tl,) +X 3 (T" kJ 4 k -TI3) +Y,E 0lJ . k 



21. Comer 6 

-2F0T;_ 1Jk + (l+2Fo)T* k = T" Jk + U 3 (t" j _, k ~T" Jk ) 
+ U 4 (T" Jk „ -T" jk ) +V 2 (T. 2 -T; k ) +V 3 (T. 3 -T" J k ) 

+ V 5 ( T -5 - T m.„) + W 2 + W 3 + W 5 + X 2 (T^ k - T 4 _ 2 ) 

+ X 3 « k -Ti 3 ) + X 5 (T:; k -T%) + Y,E 0 , Jk 



22. Perimeter 6 

(l+2F0)T* ik -2FoT;„ Jk =T" jk + U 2 (T" Jk _, + T" Jk+1 - 2T k] k ) 
+ U,(T^. k -T". k ) + V 2 (T. 2 -T” k ) + V 6 (T- 6 -T” k ) 

+ w 2 + W 6 + X 2 (T" tJ 4 k -Tt,) +X e (T; 4 t -T 4 . 6 ) + Y,E 0tJ t 
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23. Right Face 



-FoTt,.^ +(l+2Fo)T* Jk - FoT* 1JA =T"j k 
+ U 2 (T; i k ., +T" ukkl - 2T kj k ) +U 3 (T^, k -T” k ) 

+ V 2 (T. 2 -T“. k ) + W 2 + X 2 (T^ k -Tt,) + Y.Eou., 



24. Perimeter 5 



-2FoTr. I J k + (1 + 2Fo)T kJ k =T" kk +U 2 (T" J k . 1+ T" J kll -2T 
+ U 3 (t: j . 1 , -T; k ) + V 2 (T. 2 -T”. k ) +V 5 (T. 5 -T; jk ) 

+ W 2 + W 6 + X 2 (T”‘ k -Ti,) +X 5 (T^ k -Tl 5 ) + Y,E 0|Jk 



25. Corner 7 



(l + 2Fo)T- Jk -2FoT* +1 Jk "T^ + U^T^-T”.,) 
+ U.( T ",. k -, -T” lk ) +V 2 (T„ 2 -T^ k ) + V 4 (T. 4 -T" J k ) 
+ V 6 (T. 6 -T" )k ) + W 2 +W 4 + W 6 +X 2 (T^ k -Tl 2 ) 

+ X 4 (T k ; k -Tl.) +X 6 (T" 4 k -Tl 6 ) + Y,E 0 , Jk 



n 

i.j.k 
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26. Perimeter 10 



-F° T r-,j. k + «+2Fo)T* k -FoTJ„ )k =t: 0 , 

+ U s ( T I,-u - T «. k ) +u ,( T «. k -, - T „. k ) 

+ V 2 (T^ -T; j>k ) + V 4 (T- 4 -T" J k ) + W 2 + W 4 
+ X 2 (< k -Tl,) + X 4 (T^ k -Tl 4 ) + Y,E GtJ k 



27. Corner 8 

- 2FoT r-,.,. k + (1 + 2Fo)T kJ k = T kJ k + U 3 (t" j_, k - T kJ k ) 
+ U.( T r. J . k - 1 -T".0 +V 2 (T^ -T”„. k ) +V 4 (T. 4 -T” Jk ) 
+ v 5 ( t . 5 - t :,,) + w 2+ w 4+ w 5+ x 2 (t:;,- t i 2 ) 

+ X,(T" J 4 k -T; 4 ) +X s (T" J 4 k -Tl s ) + Y,E 0kJ . k 



B. Y-DIRECTION EQUATIONS 
1. Corner 1 

(1 + 2FoR ,)T f * * -2F0R,T k *; i k = T kJ k + U 4 (T kJ k+ , - T kJ k 

+ U 5 ( T ».j. k - T ,*,. k ) + V . (T-, -T" t , k ) + V 3 (T- 3 - T" k ) 

+ V 6 (T. 6 -T” J .k) + W 1+ W 3+ W 6+ X 1 (T; j 4 k -Ti l ) 

+ X 3 « k -T 4 . 3 ) +X 6 (Ty* k -Ti 6 ) +Y 1 E 0lJ k 
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2. Perimeter 3 



-FoR.T* k + (1 + 2FoR ,)T * * k - FoR.T”^ = 
t;,, + U 4 (T;, k ., - T" Jk ) + U 5 (T* +1 Jk -T* k ) 

+ V 3 (T^ -T : j k ) +V 6 (T. 6 -T" J t ) +w 3 + w 6 
+ X 3 (T"'' k -Ti 3 ) + X 6 (T" 4 k -Ti 6 ) + Y,E 0lJ . k 



3. Corner 5 



-2F0R,T;;_ l k + (l+2FoR,)T* J * k = T^ k + U .(T^ k „ -T^ k ) 

+ U 5( T r„. J . k - T *,. k ) +V 2 ( T -2- T M.k) +V 3( T -3- T “ J .k) 

+ v e(T - 6 -Ty. k ) +W 2 + W 3 +W 6 + X 2 (T" 4 k -Tl,) 

+ X 3«k - A) + X 6«k - T -e) + Y J E G , J k 

4. Perimeter 8 

(1 + 2FoR ,)T - 2FoR ,T *.,.k = T”, k + U 2 (t" , k „ + T° - 2T^ , 

+ U 5 ( T r.,j.k - T *,.k )+ v ,(T-, -T^ Jk ) +V,(T. e -T" Jk ) 

+ W 1 + W 6 + X,(T: j 4 k -T 4 .,) + X 6 (T; j 4 k -Ti 6 ) + Y 1 E 0 , Jk 
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5. Back Face 



- FoR . T “-,.k + (1+2 FoR,)T” - FoR ,T* * h - T”. k 

+ U 2 (T: j . l . 1 +T; jktl -2T; jk )+U 5 (T*, 1Jk -T* k ) + 

+ v s (T-e - T; ja ) + W 6 + X 6 (t”\ - t! 6 ) + Y, E 0 tJ k 

6. Perimeter 6 

-2 FOR ,T + (1+2FoR,)T" - T^ + U^T^.^^.-ZT^ 
+ ^5 (T T+i.j.k -T*. k )+V 2 (T. 2 -T; )k )+V 6 (T. 6 -T” Jk ) 

+ W 2 + W. + X^T”^ -Tl 2 ) +X 6 (T" ) 4 k -Tl 6 ) +¥,£ 0 ^ 

7. Comer 3 

(l + 2FoR,)T* k - 2FoR jT* j+1 k =T" J k + U 4 (t" j k _, -T" Jk ) 

+ U s (T* 1J . k - T*. k )+V,(T. 1 -T; jl ) + V 4 (T. 4 -T" Jk ) 

+ V s( T -e - T M.k) + W, + W 4 + W 6 +X,(T:* k -Tl,) 

+ X Yt" 4 . -t! 4 ) +X K (T n ' 1 . -Ti_) +y.e c ... 

4\ i.j.k 4/ 6\ l.j.k 6/ 1 ° l.j.k 
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8. Perimeter 1 



- F ° R , T I/-,.,, + (1+2FoR,)T” - FOR.T** t = T," Jk 

+ U 4 (T 1 n Jk _ J — T ” J>k ) + U 5 (T t - T * k ) + V, (T. 4 - T" „) 

+ v e(T- 6 -T J,,) +W 4 +W 6 + X 4 (T" J 4 k -Tl 4 ) 

+ X e«u- T 4 - e )+ Y , E o lJ . k 

9. Comer 7 

-2FOR.T** , + (1 + 2FoR ,)T *♦„ = k + U 4 (t", k ) 

+ U 5 ( T r.,,. k - T ', k )+ V 2 ( T - 2 - T M. k )+ V 4 (T- 4 -T" )k ) 

+ V„ (T.„ - T” k ) + W 2 + W 4 + W„ + X 2 (T” 4 k - Tl 2 ) 

+ X ,( T ," J 1 k - A) + X «« k -T 4 .,) + Y ,Eo tJ . k 

10. Perimeter 12 

(1 + 2FoR j)T * „ -2FoR,T*. I k = T" k k + U 4 (T^ k>! - t”, k ) 

+ F0(T* 1J k +T* 1Jk - 2T* k ) +V,(T., -T” k ) + V 3 (T_ 3 -T" t 
+ W, + W 3 + X 1 (T;‘ k -Ti,) +X 3 (T" 4 „ -Ti 3 ) + Y,E 0|Jk 
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11. Bottom Face 



- FoR , T tA.k + (1 + 2FOR ,)T * * - FoR.T* =T” 



J.k 



+ U 4 (T”. ktl -T” k )+F0(T*., Jk + T*, Jk -2T* k ) 
+ V 3 (T -3 - T”y k ) + w 3 + - Ti 3 ) + Y,Eo, Jk 



12. Perimeter 9 



-2FoR,T” k + (1+2F0R,)T“ = T^ + U^.-T^) 

+ Fo(T,* 1 Jk +Tf tlJk -2T* k )+V 2 (T. 2 -T: jt ) + V 3 (T. 3 -T^ k ) 
+ W 2 +W 3 + X J (T" k ‘ k -Tl 2 )+X 3 (T^ k -T 4 . 3 ) +Y,E GlJ . k 



13. Left Face 



(1 + 2FoR,)T; j * - 2FOR.T*;,., = T^ k + - 2^, 

+ Fo ( T *,. ) .k +T - 2T* k )+V,(T., - T kJ k ) 

+ W 1+ X 1 (T” 4 k -Tl,)+Y 1 E 0lJ>k 



14. Interior 



- FoR , T :/-,. k + (l+2Fo R ,)T 1 * J * k - FoR.T*,, „ =T" J t 
+ U 2 ( T M.k-i + T:. J . k .,-2T”. k )+Fo(T* 1J . k +T* 1Jk -2T*. k ) 



+ Y l E <M.J.k 
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15. Right Face 



-2F0R.T*:,,, + (1 + 2FoR ,)T ** k = + U 2 (t" ma ., + - 2T" tJ k 

+ Fo ( T *-,.ja +T r„.u -2T‘, k )+V 2 (T. 2 -T; jk ) 

+ W 2 + X 2 (T" J 4 k -Tl 2 ) + y,Eo 1J . k 

16. Perimeter 11 



(l + 2FoR 1 )T** -2FoR,T*; i k = T kj k + U 4 (T kJ k _, - T kJ k ) 
+ F °(T*,., k +T r., J .k - 2T*. k ) +V,(T., -T; jk ) 

+ V 1 (T. 4 -T;j k )+W , + w 4 + X,(T” 4 k -Ti 1 ) 

+ X 4( T r.i' k - T -4) +Y 1 E o, J . k 



17. Top Face 



— FoR jT *,*j k + (1 + 2FoR ,)T ** - FoR.T* k =T", 



j. k 



+ U 4 (T” J , k .,-T" Jk ) + Fo(T* 1Jk + T -2T* k ) 
+ V 4 (T_ 4 -T“. k ) + W 4 +X 4 (t; i \ -Ti 4 ) +Y,E 0iJ . k 
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18. Perimeter 10 



-2F0R.T*!,.,, + (1 + 2FoR ,)T ** k = T^ k + U 4 (T^., - 
+ F0(T* 1J , t + T * +1 , k “ 2T* jk ) + V 2 (T-2 - T" jk ) 

+ V 4 (T. 4 -T: jk ) + W 2 + W 4 + X 2 (T" y 4 k -Tl 2 ) 

+ X ,( T ” J 4 k - T -.) +Y , E o, Jk 



19. Corner 2 



(1 + 2FoR j)T *» k - 2FoR ,T + U 4 (t" u k41 - 
+ U 5 ( T *,.j. k - T* „ ) +V,(T., -T:, k ) +V 3 (T. 3 -T" Jk ) 
+ V 5 ( T - s - T y.,.) +W 1 +W 3 +W 6 + X 1 (T" 4 k -T 4 . 1 ) 

+ X 3 « k - t! 3 ) +X 5 (T: 4 k -Tl 5 ) + Y,E 0lJk 




20. Perimeter 4 

-For .T'j*, k + (l+2FoR 1 )T* J * k - FoR,T,*; i k = T" k 
+ U .( T ",. k ., -T". k ) +U 5 (T-,j k -T* k ) +V 3 (T. 3 -T“ k ) 
+ V 5 (T. s -T” Jk ) + W 3 + W 5 +X 3 (T” 4 k -T 4 3 ) 

+ X 5 ( T M 4 k - T - s ) +Y . E o M . k 
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21. Corner 6 



-2FoR,T” l k + (1 + 2FoR ,)T ** k = T" Jk + U^T", k „ -T", k ) 

+ U s (T*, Jk - T ,*. k ) +V 2( T - 2 - T M. k ) +V 3( T - 3 - T ",.k) 

+ V s ( T -b - T w. k ) + W 2 + W 3 + W s + X^T^ - T 4 ,) 

+ X 3 « k - T 4 . 3 ) + X s (T” 4 k - Tl 5 ) + Y ,E 0)Jk 

22. Perimeter 7 

(1 + 2FoR ,)T** k - 2FoR jT^*, k = T” u k + U 2 (T^ k ., + k „ - 2T", 
+ U 5 (T*,., k - T* k ) + V, (T-, -T” Jk ) +V 5 (T. 5 -T" Jk ) +W 1 + W 5 
+ X.« k -T 4 -,) +X 5 (T" 4 „ -Ti s ) +Y 1 E 0tJ t 



23. Front Face 



- FoR j;A. k - ( l+2FoR ,) T “ k - FoR , T r.;,. k = T r.j.. 

+ ^ 2 + T"j, k *, - 2T kJk ) + U S (T * 1Jk -T* k ) 

+ V b ( T - 5 - T m.k) +W 5 +X 5 «u-Tl 5 )+AEo tJ . k 
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24. Perimeter 5 



-2FoR,T“. k + (1 +2FoR,)T** k = U 2 (T" + T^, - 2T” , 

+ U 5 (T* 1j .„ -T* k ) + V 2 (T_ 2 -T" Jk ) +V s (T. 5 -T" Jk ) +W 2+ W S 
+ X l « k - Ti 2 ) +X 5 « k -Tl 5 ) + Y, E g ,., k 



25. Corner 4 

(1 + 2FoR j)T* * k -2F0R,T- k = T kJ k + U 4 (f g - T" k ) 
+ U ! (T* 1J . k - T* k )+V,(T. 1 -T" lJk ) +V 4 (T. 4 -T” Jk ) 

+ V (T. 5 - Ty k ) + W, + W 4 + W 5 + X,(T^ 4 t - Tl,) 

+ X 4 (T”‘ k - Tl 4 ) + X s (T”J| k - Tl 5 ) + Y,E 0 k , k 

26. Perimeter 2 

- F ° R , T ”,. k + (l + 2FoR,)T* J * k - FoRjTy* k = T" Jk 
+ U 4 (T” J . k . 1 -T," Jk ) + U 5 (T* 1Jk -T* k ) + V 4 (T. 4 -T n lJk ) 

+ V =( T -5- T "j. k ) +W 4 +W s + X .( T " J '! k -T; i ) 

+ X s« k -Ti 5 ) + Y 1 E 0lJ , 
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27. Corner 8 



-2FoR,T“ l k + (H^FoRjT," = „ + U 4 (T^ k ., - k 

+ U 5( T *, J . l - T ,*,.k ) +V 2 ( T -2- T ^) +V .( T -4- T " j .„) 

+ V 5 (T_ s -T" J k ) + W 2 + W 4 + W 5 + X a (T°\ -Tl 2 ) 
+ x 4 (Tr; k -T! 4 ) + x 5 ( T :; k - T4 . 5 ) + Y ,E 0lJk 

C. Z-DIRECTION EQUATIONS 

1. Comer 1 

(l + 2FoR 2 )T” , k -2FoR 2 T"; , k „=T** + U 3 (T" t - T,” 

+ U 5 (Tr. 1 ., k -T* t )+V 1 (T.-T“ )+V 3 (T. 3 -Tt J * k ) 

+ V a(T- 6 -T; j * k )+W 1 + W 3 + W e + X 1 (T “ 4 -T 4 .,) 

+ X 3( T t J * 1 « 4 -Tt»)+X,(T*i -T 4 e ) + Y,Ec ( j k 

2. Perimeter 8 

-FOR XIl-, + (1 + 2 FOR 2 )T" , k - FOR a T^ = T* * 

+ U 3 (T* :,. k - T* * ) + U 5 (T *, 0 k - T * k ) + V,(T., - T” 

+v 6 (t. 6 -t“ ) +w, + w t +x,(T**; -Tl.) 

+ X s ( T " k 4 -T 4 - 6 )+Y 1 Eo,. J>k 
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3. Comer 3 



-2FoR a T" , k . 1 + (l+2FoR s )T"J = T 1 V k + U 3 (T** k - T**„ ) 

+ U 5 (Tr.,, k - T*. k ) +V 1 (T- 1 -T“ ) +V 4 (T. 4 -T“ ) 

+ V 6 (T. 6 -T“ ) + W, + W 4 + W 6 + Xj(T*j* k -Tl.) 

+ X .( T M* k 4 - T4 -«) + X 6( T ,** k -T^+Y.Eo,.,,, 

4. Perimeter 12 

(1 + 2FoR 2 )T“* k - 2FoR 2 T"* k+I = T* k + U 3 (T 1 * ]4 * lk - T" ) 

+ Fo ( T ,-..,. k +T r..j. k - 2T *,. k ) + V 1 (T.-T** )+V 3 (T. 3 -T” ) 

+ W, +W 3 + X,(tJ j * 4 -Ti 1 ) + X 3 (T 1 ^ ”T- 3 ) +Y,EG 1Jk 

5. Left Face 

-2 FoR 2 T"j k _ I + (1 + 2 FoR 2 ) T "* k = T,* k + U 3 (T**, k - T" ) 

+ F°(T*->*k +T*,, k - 2T* k ) +V,(T- 1 -T” ) + V 4 (T. 4 -T** k ) 

+ W,+W 4 +X 1 (T** k -Ti,) + X 4 (T** k -T; 4 ) + Y,Eo 1Jk 
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6. Perimeter 1 1 



-2FoR ! T-.,+(l + 2FoR 2 )T-=T** k + U 3 (T* * k - T** ) 

+ Fo ( T r-,. J . t + T*,.,. k - 2T*, k ) +V,(T., -T” ) + V 4 (T„„ -T** ) 

+ Wj +W 4 +Xj(T** k _ t1j) + X 4 (T*j* k -Tl 4 )+y,Eo 1 . J . k 



7. Comer 2 



(1 + 2FoR 2 )T"* k - FoR 2 T"’ k+] = T ( * * + U s (T**. k - T>* ) 

+ U 5 ( T *,,. k -T*. k )+V,(T.,-T** k )+V 3 (T. 3 -T” ) 

+ V 5( T -5- T “ k ^W. + W^W^X^T” 4 -T-,) 

+ X 3 (T“ 4 -Ti 3 ) + X 5 (T-; k 4 -Tl 5 ) + YjEg, Jk 



8. Perimeter 7 

- F 0 R 2 T" k ., +(1 +2F0R 3 )T"' k - FoR 3 T;.;' k ., = T** k 
+ u.fr;:,, -T** k )+o 5 (T 1 *. l . ) . k -T t * k ) 

fV ,(T-,-T“ )+V 5 (T,-T>* )+W 1 + W 5 

+x ,( T :; k 4 - Ti -,) +x s ( T “ k -T: 5 )+y,Ec,., k 
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9. Comer 4 



-2 FoR 2 T"* k _ I + (1 + 2FoR 2 )T"‘ = T“ 


+ u 3 (t** k 


+ U 5 (T* 10 , k - T* „ 


)+V,(T- 1 -T“ 


)+v 4 (T. < -T 


+ V s(T-5 _T t)* k ) 


tWj + W. + WjtX, 


/ 'p * * 4 T 4 ^ 

l v i.j.k 1 “ 1 ) 


+ X 4( T " k -T-.) 


+x.(T-{ -Tl,) 


+ Y i Ec i.j.k 



10. Perimeter 3 

(1 + 2FOR 2 )T";' k - 2 FOR .T”; 1 ,,, - T ♦ * 

+ U,(T** k +T “,. k - 2T,V k ) + U 5 (T* 1Jk -T* >k ) 
+ V s(T-3-Tt; k )+V 6 (T- 6 -T” ) + W 3 + W 6 

+ X 3( T " k - T - 3 ) + X «( T " k -Tl.J+Y.Ec.j,, 

11. Back Face 

-FoR.T";!-. + (1 + 2FoR 2 )T" k ‘ k - FoR 2 T"J' k41 =T; 

+ U,(T**. k + T*;, k - 2T** ) + U s (T* uk -T* k ) 

+ V 6 (T- 6 -T“ )+w 6 +x 6 (t;*‘ -O+Y.Eoy., 



.J.k 



.n+1 



. n+1 



— *T * * 
i.j.k 
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12. Perimeter 1 



-2FoR 2 T"*' k _ l + (1 + 2 FoR 2 )T”: , |1 = T*< 

+ U,(T”. k +T** k - 2T** k ) + U 5 (T* |Jk -T* k ) 

+ v 4 (t. 4 -t: j * )+v 6 (T. 6 -T;; k ) + w 4+ w 6 

- T ‘-4)+X 6 (T* J *; - T- 6 ) + Y,EG 1Jk 

13. Bottom Face 

(1 + 2FoR 2 )T"* k - 2FoR 2 T"' k4 , =T* * 

+u ,( T t;,. k + T r;,. k - 2T ,” ) 

+ Fo (T,“. k - 2T** k )+V 3 (T_ 3 -T; j * k ) 

+ W 3+ X 3( T r.,* k - T4 -3) + Y , E <=,,. k 



14. Interior 

+ (1 + 2FoR 2 )T"* k - FoR 2 T"' k41 = T" 
+ U 1 (T“ k + T* J4l k - 2T** k ) 

+ Fo(T l * 1)k + T* +I)k - 2T* k ) + V 4 (T. 4 -T“ ) 

+ W 4 + X 4 (T * * k -Tl 4 ) + Y,Eo, J k 
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15. Top Face 



-2FoR 2 T"‘ k -, +(! +2FoR 2 )T" , k = T** 

+ u ,(T“,. k + T** „ - 2T-- ) + Fo(T*, Jk +T* 1)k -2T* k ) 

+ v 4 (t. 4 -t; j * )+w 4 + x 4 (T;; k 4 -t 1 4 ) +Y.EC,.,., 



16. Perimeter 4 

(1 + 2FoR 2 )T"* k - 2FoR 2 T"‘ k41 =T* * 

+ U,(T“,. k +T“„ - 2T” ) + U 5 (T*, Jk -T* k ) 

+ V 3 (T- 3 -T“ ) +V 5 (T_ 5 -TJj\ )+w 3+ w 5 

+ X 3 ( T “ k -Tl^ + Y.Eo,^ 



17. Front Face 

- F ° R J "L, + (1 + 2 For 2 )T"* k - For 2 T 4l - T ,** 
+ u .( T t,:,. k + T“. k - 2T" ) 

+ U s (T*, Jk -T*„ )+v 5 (t. 5 -t; j * ) 

+ W 5 +X 5 ( T r.j* k -TiO + Y.Ec,.,., 
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18. Perimeter 2 



-2FoR 2 T,J kl + (1 + 2 FoR 2 )T”j k = T,* k 
+ U,(T“ k +T** k - 2T** k ) + U 5 (T* 1Jk -T* k ) 

+ V 4 (T-„-T** k ) +V 5 (T. 5 -T** k )+W 4+ W 5 

+x ,( T " k 4 - T4 -4) +x 5 ( T “ k -TlJ+Y.Eo,^ 

19. Comer 5 



(1 + 2FoR 2 )T" J k - 2FoR 2 T, j k41 = T * k 
+ U 3 (T”,. k - T** )+U 5 (T* 1Jk -T‘. k ) 

+ V 2 (T. 2 -T** k )+V 3 (T^-Tt; k )+V 6 (T. 6 -T“ ) 

+ W 2 + W 3 + W 6 + x 2 (t ** k -t! 2 ) 

+ X.,(T” k 4 -Tl 3 )+X 6 (T* J * k -Tl e ) +Y,Eo 1)k 

20. Perimeter 6 

- FoR 2 T ”.U +0 + 2F0R 2 )T;;‘ k -FoR a T“‘ k41 =T*; k 

+ U 3 (T;A. k - T ” k ) +U .( T * 1J . k - T** k ) 

+ V 2 (T. 2 -T“ )+V 6 (T.,-T** )+w 2+ w 6 

+ x 2 (T” 4 -Tl 2 )+X 6 (T* J *‘ ~t 1 6 ) +YjEG IJk 
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21. Comer 7 



-2FoRj”i_ 1 + (l+2FoR 2 )T^ = T** + U s (T*;,. k - T** ) 



+ U 5 ( T - T ij. k 


) + V 2 (T- 2 -T“ 


) + V 4 (T. 4 -T" 


+ v 6 (t. 6 -t; ] * ) 


+ w 2 + w 4 +w 6 +x 


2 ( T ,V k - T - k ) 


+ x 4 (T ” k 4 -t 4 - 4 ) 


+x 6(T,*; k -A) 


+ Y l EG i.J.k 



22. Perimeter 9 

(1 + 2FoR 2 )T”;' k - 2FoR 2 T"‘ ll ., =T“ 

+u s (TrA. k - T “„ ) +Fo ( T r.,. J . k +T ,:,,. k - 2 T r,. k ) 

+ V a (T-j-T“ ) +V 3 (T . 3 -T“ ) + W 2 + W 3 

+ x 2 (t “„ 4 -Tl 2 )+X 3 (T“ k -Ti 3 ) +Y 1 Eo 1Jt 

23. Right Face 

-F0R 2 T-., + ( 1+ 2F0R 2 )T- -FoR 2 T"'„ =T" 

+ U 3 ( T M*.. k - T " k ) +F °( T *, J . k +T* 1Jk -2T* t ) 
+ v 2 (T- 2 -t; j * )+w 2 + x 2 (t“ 4 -Tl 2 ) + Y,Eo, J . k 
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24. Perimeter 10 



-2FoR 2 T";L, + (1 +2 FoR 2 )T^ = T** + U,(T - T** ) 

+ Fo(T*, Jk +T r.,j. k -2T* k )+V 2 (T- 2 -T“ ) 

+ V 4 (T.,-T” )+W 2 + w 4 +x 2 (T ; J * k 4 -t 4 . 2 ) 

+ X 4 (T *; k 4 -T 4 . 4 )+Y 1 Eo tJ . k 



25. Comer 6 

(1 + 2FoR 2 )T"‘ k - 2 FoR 2 T;.; , ,„ 1 = T» * 

+ U 3( T ”i.k - T** )+U 5 (T‘ IJk - T* k ) 

+ v 2 (t. 2 -t” ) +v 3 (t. 3 -T“ )+V 5 (T. 5 -T” ) 

+w 2+ w 3 + w 5+ x 2 (t* j - 4 -t 4 - 2 ) + x 3 (t” 4 -t 4 _ 3 ) 

+ X 5 ( T "k -T 4 . 5 )+Y,Ea, J . k 



26. Perimeter 5 

-FoR 2 T”* k _, +(1 + 2F0R 2 )T"' k -FoR 2 T” 4, k4l = T** k 

+ U s( T tA.k - T*‘ k )+U s (T* 1Jk - T” ) 

+ V 2 (T. 2 -T“ )+v 5 (t^-t” )+w 2+ w 5 

+ X 2 ( T “k -t! 2 )X 5 (t *; 4 -Tl 5 ) +YjEc lkk 
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27. Comer 8 



-2 FoR + (1 + 2 FOR 2 ) = T * * 


+ u 3 ( T tA. k - T " k 


+ U 5 (T*„„ - T,* k 


) + v 2 ( T » 2 - T *; k 


) + V.(T- 4 -T“ ) 


+ v 5 (t- 5 -t; j * ) 


+ w 2 + w 4 + w 5 + X 


( T * * 4 - T 4 ^ 

2\ 1 f.j.k 1 “2 / 


+ X ,( T “* 4 - T -0 


+x (t** 4 -t 4N ) 

5 \ J.J.k 1 - 5 / 


+ Y l EG t.j.k 
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APPENDIX C 



PROGRAM EXPLICIT 



ft ft ft. ft ft ft ft ft ft * ft * ft * * ft ft ft ft * * * * * * * ft ft * * ft * ft. ft * * $ * ft. * £ ft ft ft ft. ft. * £ <« ft ft. * * * ft ft. ft ft. ft. ft ft ft ft * ft * ft ft. ft * ft 



* PROGRAM EXPLICIT 

* 

ft ft ft * ft ft ft ft. ft * ft. * * * ft * * * * * ft * * * * ft * * * * * ft * * * * ft * ft * ft * * * * ft * * * * ft £ * * * ft * * * * * ft ft ft * ft * ft * * ft 



* EXPLICIT METHOD FOR TRANSIENT HEAT CONDUCTION IN THREE 

* DIMENSIONAL CARTESIAN COORDINATES 



* ft ft * ft * ft. * £ * ft £ * * * * ft ft ft ft * £ * * ft * 3}? * * * ft ft. ft ft ft ft ft ft ft * ft * * * * $ * * * * ft * ft $ * * ft. ft. ft ft. ft ft ft ft ft ft ft * * * * 

* * ft ft ft * ft ft * * * * ft * A * * $ * * ft £ * ft * ft $ if- * ft ft * * ft ft * * * * * ft A * A * * ft * ft. ft ft. ft ft ft. ft. ft- ft ft ft. ft. ft. ft. ft. ft. ft ft. ft ft. ft ft ft. 

ft 

* DESCRIPTION OF VARIABLES USED IN PROGRAM EXPLICIT 



* NOTE: SI UNITS ARE USED IN THIS PROGRAM UNLESS 

* OTHERWISE NOTED 

ft 

* 1. USER SPECIFIED VARIABLES 

ft 

* A. GEOMETRIC DIMENSIONS: 

ft 

* DELX-DI STANCE INCREMENT IN THE X DIRECTION 

* DELY-DI STANCE INCREMENT IN THE Y DIRECTION 

* DELZ-D I STANCE INCREMENT IN THE Z DIRECTION 

* 

* I-— X COORDINATE OF NODE 

* J— -Y COORDINATE OF NODE 

* K— -Z COORDINATE OF NODE 



LX— LENGTH OF MODEL IN THE X DIRECTION 
LY— LENGTH OF MODEL IN THE Y DIRECTION 
LZ— LENGTH OF MODEL IN THE Z DIRECTION 

M— NUMBER OF NODES IN THE X DIRECTION 
N — NUMBER OF NODES IN THE Y DIRECTION 
P — NUMBER OF NODES IN THE Z DIRECTION 



B. MATERIAL PROPERTIES: 

CP— SPECIFIC HEAT 

Kl— THERMAL CONDUCTIVITY 

RHO-DENSITY 



C. INITIAL DATA: 



FLUX-SURFACE FLUX 

FREQ-NUMBER OF TIME STEPS BETWEEN SUCCESSIVE 
TEMPERATURE PRINTOUTS 

II HEAT TRANSFER COEFFICIENT 

TAMB-AMBIENT TEMPERATURE 
TINTT-INITIAL TEMPERATURE OF ALL NODES 
VAL— LOGIC NUMBER TO DETERMINE IF VALIDATION 
TEST IS TO BE RUN 



2. NON-USER VARIABLES 



A. GENERAL: 



BIO— BIOT NUMBER 

COUNT-COUNTER FOR DETERMINING THE NUMBER OF ELAPSED 
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TIME STEPS 



DELT— TIME INCREMENT 

FO FOURIER NUMBER 

TIME— DURATION OF CALCULATIONS 

TNEW— MATRIX OF TEMPERATURES AFTER TIME STEP DELT 



* B. EXTERNAL PROGRAMS: 

❖ 

* CTIME USED IN EXTERNAL IMSL ROUTINE FOR DETERMINING 

* AMOUNT OF CPU TIME USED 

* N'OU'T USED IN EXTERNAL IMSL ROUTINE FOR DETERMINING 

* AMOUNT OF CPU TIME USED 

* TIMEO CPU TIME AT BEGINNING OF CALCU LATIONS 

* TIM El CPU TIME AT END OF CALCULATIONS 

* UMACH USED IN EXTERNAL IMSL ROUTINE FOR DETERMINING 

* AMOUNT OF CPU TIME USED 



# * * * * * * * * * * * 'I: # * * * * -V * * # A £ Ss * * * * * . J .: * * * * * * * * * * .*5: * * # * * # * * # * * * * * * * * * * * 

* * * * * * * * * * * * * * * * * * * * * * * * * * -'I'- * * * * * * A * .'1: * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 



DIMENSION T(25.25.25). TNE\V(25.25.25) 



REAL TAMB. H. RIIO. ALFA, FO. FLUX, R. Rl, R2, K1 

REAL DELX. DLLY. DELZ, DELT. TIME, TINIT, TNEW, BIO 

REAL TDIFE1, TDIFE2, TDIFE3 

REAL LX, LY, LZ 

REAL CTIME, TIMEO, TIME1 

INTEGER I, J, K, M, X, P, COUNT, FREQ, VAL. NOUT 
EXTERNAL CTIME, UMACH 

COMMON TNEW, TIME. DELT, M, X, P, TDIFE1. TDIFE2, TDIFE3. TINIT 
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COMMON COUNT, FREQ 



* 

* 4 * * * * 4 4 * * 4 4 4 4 4 4 4 * 4 * * * 4 * 4 * * * 4 * * 4 4 4 * * 4 * 4 * * * * * * * * * * 4 * * * * * * * * * * * * * * * * * 

* * ■:■■ * * * * * * * * 4 * * * * * « * * * * * * * * * ■•:•■ * * * * 4 * * * * * * 4 * ■:• * * * * * * * •:•■ 4 * 4 * * 4 * * .4 * 4 * * 4 4 4 * * 4 * 
4 

* INITIALIZING CPU TIME 

* 

TIMEO = CTIMEO 

4 

4 4*4* 4 4 4 4 4 4 * 4 * 4 4 4 4 4 4 4 4 4 4444 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 44444444 4 * 4444 4 4 4 44*44 4 4 4 
4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 * 4 * 4 4 4 4 4 4 * 4 4 4 4 4 4 4 4 44 4 4 4 4 4 4 4 4 4 4 * 4 4 * 4 4 4 4 * 4 4 4 * 4 4 4 4 4 4 4 4 4 4 4 



* VARIOUS PARAMETERS TO BE USED THROUGHOUT THE PROGRAM 

* WILL NOW BE READ INTO THE PROGRAM. 



WRITEC*,*) 'ENTER TAMB. THE AMBIENT TEMPERATURE.' 

REA DC®.*) TAMB 

WRITE!*,*) ENTER TINTT. THE INITIAL TEMPERATURE.' 
READ(V) TINTT 

WRITE(*.*) 'ENTER IT. THE HEAT TRANSFER COEFFICIENT.' 
READ(*,*) II 

WRITE!*.*) 'ENTER Kl. THE THERMAL CONDUCTIVITY.' 
READ(*.*) Kl 

WRITEC*.*) ENTER CP. THE SPECIFIC HEAT OF THE MATERIAL.' 
READ(*.*) CP 

WRITEC*,*) ENTER RIIO, THE MATERIAL DENSITY.' 

READ(*,*) RHO 
WRITE)*,*) 'ENTER FLUX.' 

READ(*,*) FLUX 

WRITEC*,*) ENTER LX, THE LENGTH IN THE X DIRECTION.' 
READ(v) LX 

WRITEC*,*) 'ENTER LY, THE LENGTH IN THE Y DIRECTION.' 
READC*,*) LY 

WRITEC*,*) ENTER LZ. THE LENGTH IN' THE Z DIRECTION.' 
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READ!*,*) LZ 

WRITE!*.*) 'ENTER M, THE NUMBER OF NODES IN THE X DIRECTION.' 
READ!*,*) M 

WRITE!*,*) 'ENTER N, THE NUMBER OF NODES IN THE V DIRECTION.' 
READ(*,*) N 

WRITE!*,*) 'ENTER P, THE NUMBER OF NODES IN' THE Z DIRECTION.' 
READ!*,*) P 

WRITE!*,*) 'ENTER FREQ, THE NUMBER OF TIME STEPS BETWEEN" 
WRITE!*,*) 'SUCCESSIVE PRINTINGS.' 

READ!*,*) FREQ 

WRITE!*,*) 'ENTER VAL: IF YOU WANT THE VALIDATION SUBROUTINE' 
WRITE!*.*) CALLED, USE 1; IF THE VALIDATION SUBROUTINE IS NOT' 
WRITE!*,*) DESIRED. USE 0.' 

READ!*.*) VAL 



* * s-. * * * ft * * * * * ft ft ft * ft ft ft ft * ft ft ft * ft ❖ ft * * ft ft ft ft ft ft ft •> ft ft ft ft ft * ft * * * ft * * ft ft ft ft >:•. ft * ft « * ft ft ft ft ft * ❖ * 

* * ft * ft * * * * * * ft * ft ft * * * ft ft ft ft * ft * * * * * * ft s-. * * ft ft « « ft * * * * * ft * * * ft * * ft * * ft ft ft ft ft ft ft * « * * * * * 



* CALCULATIONS OF CONSTANTS TO BE USED THROUGHOUT THE 

* PROGRAM 



COUNT = 0 
TIME = 0.0 



DELX = LX (M-l) 

DELY= LY (N-l) 
DELZ=LZ (P-1) 

R= DELX (DELY*DELZ) 
R1 = (DELX**2) (DELY**2) 
R2 = <DELX**2) (DELZ**2) 

ALFA= KL(RHO*CP) 

B I O = II Kl -DELX 
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F0= 0.5 (1 + R1 + R2 + BIO*SQRT(R2)) 

DELT = FO*(DELX°*2) ALFA 

* NOTE: THE FOLLOWING EXPRESSION FOR DELT WAS USED FOR 

* VALIDATING THE EXPLICIT TECHNIQUE WITH THE CLOSED-FORM 

* SOLUTION. 

* 

* DELT = 0.125*(DELY**2)'ALFA 

ft 

* * * ft ft ft * * * * * * * ft. ft. * * * * * * * * * * * * * * * * * * * >:-. ft ft * ft * * * ft * * * * * * * * * * ft ft ft ft « ft * « * ft « * ft ft <•. ft 

« « ft ft « ft ft « ft « ft « ft ft * ft ft « « « « ft * ft * ft * ft * « * * * * * * ft ft. * * ft. « * « « ft. * « « * « « ft ft. « ft ft « « ft. ft. <•. ft. * * * * * 

ft 

* INITIALIZE ALL TEMPERATURES TO TINIT, THE INITIAL 

* TEMPERATURE. 

* 

DO 10 1= 1. M 
DO 15 J= 1. N 
DO 20 K= 1. P 

TNEW(IJ.K)= TINIT 

20 CONTINUE 
15 CONTINUE 
10 CONTINUE 

* 

* 

* * * * * * * g * * * * * g * * * g * * * * * * * * * * g * * * g * * * * g g g * * * g * * * * * * * * * g * * * * * * * * * * * * * * * * * 

* * g * * * * * * * g * * * g * g g g g g g g g g g g * gggggggggggggggggggggggggggggg g ggggggggggg g 

* 

* THE OLD TEMPERATURES WILL NOW BE SAVED FOR FUTURE USE 

* IN THE TEMPERATURE EQUATIONS. 

g 

DO 30 1= 1, M 
DO 35 .1=1, N 
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DO 40 K = 1. P 



T(I.J.K) = TNEW(I.J.K) 

40 CONTINUE 
35 CONTINUE 
30 CONTINUE 

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * >> >> * * * * * * * * * * * * * * * * * * * * ****** 
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ************** * * * * * * 

INCREMENTING TIME AND COUNTER 

1 COUNT = COUNT + 1 
TIME = TIME + DELT 

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

NODE EQUATION'S TOR THE THREE DIMENSIONAL BLOCK ARE 
NOW PRESENTED. 

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

INTERIOR NODES 

DO 50 1 = 2, M-l 
DO 55 J = 2. N-l 
DO 60 K = 2, P-1 

TNEW(1,J,K)= T(I.J,K) + Rl*FO*( T(1,J-1,K) + T(I,J + l.K) - 2* 

& T(LJ.K) ) + R2*FO*( T(I,J,K + 1 ) + T(IJ,K-1) - 2*T(I,J,K) ) 
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& + FO*( T(I + 1.J.K) + T(I-U.K) - 2*T(IJ.K) ) 



60 CONTINUE 
55 CONTINUE 
50 CONTINUE 



* * * * * * * .•>. 4 4 * * * * * # * * * * * * * * 4 >:■ * * * 4 * « « * * * « 4 44 4 * * * * * * * * * * * * * * * * 4 * * * * * * * * * 4 * 
4 4 4 4 4 * 4 * 4 * 4 * * * * & * * 4 4 * 4 4 4 * 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4444444444444 4 4 4 4 4 4444 44 



* CORNER NODES 

4 

* CORNER NR. 1 



1 = I 
J= 1 
K= I 

TNEW(I.J.K)= T(I.J.K) + 2*R2*FO*( T(I,J.K + 1) - T(1.J,K) ) + 
& 2*Rl*FO*( T(I.J + l.K) - T(U,K) ) + 2*F0*( T(I+ 1,J.K) - 
& T(I.J.K) ) + 2 * F O * R 1 * D E L Y * F L U X K1 



- CORNER NR. 2 



I = M 
J= 1 
K = 1 

TNEW(I,J,K)= T(I,J,K) + 2*R2*FO*(T(I,J,K+l)-T(I f J f K)) + 
& 2*Rl*FO*( T(I,J + 1,K) - T(I,J,K) ) + 2*FO*( T(I-1,J,K) - 
& T(I,J,K) ) + 2*F0*R1*DELY*FLUX,'K1 
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* CORNER NR. 3 



<5 



1=1 
J= 1 
K = P 

TNEW(1.J,K) = T(1,J,K) + 2*FO*( T(1 + 1,J,K) - T(1.J,K) ) + 
& 2*Rl*FO*( T(1,J + 1,K) - T(I,J,K) ) + 2*R2*FO*( T(I,J,K-1) 
& T(1,J,K) ) + 2*FO*Rl*DELY*FLUX K1 



CORNER NR. 4 



I = M 
J= 1 
K=P 



TNEW(IJ.K) = T(U.K) + 2*FO*( T(I-1.J,K) - T(I.J.K) ) + 

& 2*Rl*FO*( T(1.J + 1.K) - T(l.J.K) ) + 2*R2*FO*( T(1.J,K-1) 
& T(l.J.K) ) + 2*FO*Rl*DELY*FLUX K1 



* CORNER NR. 5 

1=1 
J = N 
K= 1 

TNEW(1,N,1) = T(l,N.l) + 2*FO*( T(2 ,N,1) - T(1,N\1) ) + 
& 2*Rl*FO*( T(i,N-l.l) - T(1,N\1) ) + 2*R2*FO*( T(1,N, 2) 
& T(1,N,1) ) + 2*F0*B10*(R1 **0.5)*( TAMB - T(1,N,1) ) 
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* CORNER NR. 6 



M* 



I = M 
J = N 
K = 1 

TNEVV(I,J,K) = T(I,J,K) + 2*FO*( T(1-1,J,K) - T(I,J.K) ) + 

& 2*RI*FO*( T(IJ-1,K) - T(I.J.K) ) + 2*R2*FO*( T(I,J,K+ I) - 
& T(I,J,K) ) + 2 * FO * B 1 0*( R 1 * *0. 5 ) *( TAMB - T(I,J,K) ) 



* CORNER NR. 7 

1=1 
J = N 
K = P 



TNEW(I.J.K) = T(l.J.K) + 2*FO*( T(l + EJ.K) - T(I.EK) ) + 
& 2*Rl*FO*( T(l.J-l.K) - T(EJ.K) ) + 2*R2*FO*( T(I,J.K-I) - 
& T(U.K) ) + 2 * F O * B 1 0 * ( 1 * * 0 . 5 ) * ( TAMB - T(I.J,K) ) 



* CORNER NR. S 



1 = M 
J = N 
K = P 

TNEW(EJ.K) = T(I,J,K) + 2*FO*( T{I-1,J,K) - T(1,J,K) ) + 

& 2*Rl*FO*( T(I.J-1,K) - T(I,J,K) ) + 2*R2*FO*( T(I,J,K-1) - 
& T(I,J,K) ) + 2 < FO*BIO*(RI**0.5)*( TAMB - T(1,J,K) ) 



>!*. :J: t'fi <• 



* * ,v £ Jjt * A # A A ,V. £ Jj. A A * ,v A jV * # * A ,V £ A A £ ,J, * * A jJ. £ jV ,J. jJ, A A A 
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* * * * * * * * * * * * * * ❖ * * * ❖ * * * * * * * * * * * * * * * * * * ••:*. * * * * * * * * * * * * * # * ■:■ * * * * * * * •:- * * # * * * * * * 



* ADIABATIC SURFACES 

* 

* FRONT YZ SURFACE 

* 

DO 70 J = 2. N-l 
DO 75 K = 2, P-1 
I = M 

❖ 

TXEW(UJ.K) = T(I ,J,K) + 2*FO*( T(I-1,J,K) - T(I,J.K) ) + 
& Rl*rO*( T(I.J-I.K) + T(I,J+ UK) - 2*T(I.J,K) ) + 

& R2*FO*( T(IJ,K-1) + T(IJ.K + 1) - 2*T(UUK) ) 

5 /. 

75 CONTINUE 
70 CONTINUE 



* BACK YZ ADIABATIC SURFACE 



DO SO J = 2, N-l 
DO S5 K = 2. P-1 
1=1 



TNEW(I.J,K) = T(l.J.K) + 2*FO*( T(I+ 1.J.K) - T(I,J.K) ) + 

& Rl*FO*( T(I,J-UK) + T(I,J+ 1,K) - 2*T(I,J,K) ) + 

& R2*FO*( T(I,J,K-1) + T(U,K+ 1) - 2*T(U,K) ) 

* 

85 CONTINUE 
80 CONTINUE 

* 

* 

ftftftftftftftftft. ftftftftftftftftft* ft ftftftftftftftftft: ftftftftftftftftft ftft ft ft. it***************************** 
ft * * ft S'. * * S'. * ft * ft * ft * * * * * * * ft ft * * * ft ft * * ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft- ft * * ft ft. ft. ft. ft. ft ft ft ft. * ft ft ft ft ft ft ft ft ft ft 
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BOTTOM XV ADIABATIC SURFACE 



DO 90 1 = 2, M-l 
DO 95 J = 2, X-l 
K = 1 

TXEW(I,J,K) = T(1J,K) + FO*( T(I+ U,K) + T(I-1,J,K) - 
& 2*T(I.J,K) ) + Rl*FO*( T(U-1,K) + T(IJ+ I,K) - 2*T(I,J,K) ) + 
& 2*R2*FO*( T(I.J.K+ 1) - T(I,J,K) ) 

95 CONTINUE 
90 CONTINUE 



TOP XV ADIABATIC SU RFACE 



DO 100 1 = 2. M-l 
DO 105 J = 2. N-l 
K = P 

TNEW(I.J.K) = T(I,J,K) + FO*( T(l + l.J.K) + T(I-1,J,K) - 
& 2*T(I.J.K) ) + Rl*FO*( T(EJ-l.K) + T(I,J+ 1,K) - 2*T(U,K) ) + 
& 2*R2*FO*( T(UJ.K-l) - T(I.J.K) ) 

105 CONTINUE 
100 CONTINUE 



J: * :>• * * * *. * * * * * * * $ * * * * * * * * * * *. * * * * * * * * * * * * * * * * * * *: * * * ************ * * * * * * * 
5s * * * * * * * * * * * * * * * * * * * * * * * * * .* * * *. * * * *****:«:***:*: * * * * * * * * * ***** * * * * * * * * * * 



RIGHT XZ CONVECTIVE SURFACE 
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DO 110 1 = 2. M-l 
DO 115 K = 2, P-1 
J = X 

TNEW(I,J.K) = T(I,J,K) + FO*( T(I + 1,J,K) + T(I-1,J,K) - 
& 2*T(I,J,K)) + 2*R1*F0*( T(I.J-1,K) - T(I,J,K) ) + R2*FO* 
& ( T(I,J,K + 1) + T(1.J,K-1) - 2*T(I,J,K) ) + 

& 2 * F O * B I O * ( R 1 * * 0 . 5 ) *( TAMB - T(1,J,K) ) 

t 

115 CONTINUE 
110 CONTINUE 



s-- « s-- * * * * * * * * * * * * ,-f. * * * * * * * c. * * * sjt * * * * si * si >?. s< si s-. * » s< * * * * * * si * * * s> s-. * s * « * * * a s. s-. * s-. 

* St * S: St St * * St St St St St St St St St St St St St * St St St St St St St St St St St St St * St St St St St St S> St St * St St St St St St St St St St St St St St St St St St St St St * St St St St 
St 

* LEFT XZ FLUX SURFACE 

DO 120 1 = 2. M-l 
DO 125 K = 2, P-1 
J= 1 

* 

TXEW(LJ.K) = T(I.J.K) + FO*( T(I-1.J,K) + T(I+l,J.Kj- 
& 2*T(I.J,K) ) + 2*RI*FO*( T(I.J+ l.K) - T(I.J.K) ) + 

& R2*FO*( T(LJ.K-l) + T(I.J.K + 1) - 2*T(I.J,K) ) + 

& 2*FO*RI*DELY*FLUX K1 

125 CONTINUE 
120 CONTINUE 

* 

St 

* * St St « St St S< * St «t S- * Si tSt St Si s. S< * * * * St St St * * * * * * St Si Si * * * S= St St * # * s. * S- Si * * * St St St s. * * « * St * * * * St St * * * * * * 
St # St St St * <t St St St St St St St St St St S'. St St s: St St St St St St SI St St St St St St St St St St St St St St St St St St St St Si St St St St St St St St St St St St St St St St St St St St St St 

<t 
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OUTER PERIMETER SURFACES 



PERIMETER NR. I 



DO 130 J = 2, X-l 
1 = I 
K = P 

T\EW(I,J,K) = T(I.J,K) + Rl*FO*( T(I,J-1,K) + T(I,J+1,K) 
& 2*T(I.J,K) ) + 2*FO*( T(I+ U,K) - T(I,J,K) ) + 

& 2*R2*FO*( T(I,;.K-1) - T(I,J,K) ) 

130 CONTINUE 



PERIMETER NR. 2 



DO 135 J = 2. X-l 



I = M 



K=P 



TNEW(U.K) = T(IJ.K) + RI*FO*( T(U-I.K) + T(I,J + 1,K) 
& 2*T(1.J,K) ) + 2*FO*( T(I-IJ.K) - T(1J.K) ) + 

& 2*R2*FO*( T(I.J.K-l) - T(I,J.K) ) 

135 CONTINUE 



PERIMETER NR. 3 



DO 140 J = 2, X-I 
1=1 
K = 1 



9S 



TNEW(IJ.K) = T(I.J.K) + R1*F0*( T(IJ-l.K) + T(I,J + 1,K)- 
& 2*T<IJ.K) ) + 2*FO*( T(I + 1.J.K) - T(I.J,K) ) + 

& 2*R2*FO*( T(I J.K + 1) - T(I.J.K) ) 

140 CONTINUE 



* PERIMETER NR. 4 

* 

DO 145 J = 2, N-l 
I = M 
K = 1 

TNEW(IJ.K) = T(I.J.K) + Rl*FO*( T(I,J-1.K) + T(U+1.K)- 
6c 2*T(IJ.K) ) + 2*FO*( T(I-UJ.K) - T(I.J.K) ) + 

& 2*R2*FO*( T(U.K + 1) - T(l.J.K) ) 

145 CONTINUE 



PERIMETER NR. 5 



DO 150 K = 2. P-1 
I = M 
J = N 



TN'EW{I,J,K) = T(I,J,K) + R2*FO*( T(I,J,K + 1) + T(I,J,K-1)- 
& 2*T(I,J,K) ) + 2*R1*F0*( T(EJ-l.K) - T(I,J,K) ) + 

& 2*FO*( T(I-1,J,K) - T(I,J,K) ) + 2*FO*BIO*(R1**0.5)*( TAMB - 
& T(I,J,K) ) 

150 CONTINUE 
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PERIMETER NR. 6 



DO 155 K = 2, P-1 
1 = 1 
J = N 

TNEW(I,J,K) = T(I.J,K) + R2*FO*( T(1,J,K + 1) + T(I,J,K-1) - 
& 2*T(U,K)) + 2*Rl*FO*(T(I,J-l,K)-T(I,J,K)) + 

& 2*FO*( T(I + 1,J.K) - T(1,J,K) ) + 2*FO*BIO*(R1**0.5)*( TAMB 
& - T(I.J.K) ) 



155 CONTINUE 



PERIMETER NR. 7 

DO 160 K = 2. P-1 
I = M 
J= 1 



TNEW(I.J.K) = T(I.J.K) + FO*R2*( T(I.J.K+1) + T(IJ,K-1) - 
& 2*T(I,J,K) ) + 2*FO*( T(I-l.J.K) - T(I.J.K) ) + 

& 2*FO :: R 1 *( T(IJ + 1,K) - T(I.J,K) ) + 2*FO*Rl*DELY*rLUX K 1 

160 CONTINUE 



PERIMETER NR. 8 

DO 165 K= 2, P-1 
1 = 1 
J= 1 

TNEW(U.K) = T(I.J.K) + FO*R2*( T(I,J,K+1) + T(IJ.K-l) - 
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& 2*T(I,J.K) ) + 2*F0*( T(I + l.J.K) - T(I.J.K) ) + 

& 2*F0*R1*( T(I.J+ I.K) - T(I.J.K) ) + 2*F0*R1*DELY*FLUX K1 

165 CONTINUE 



* PERIMETER NR. 9 



DO 170 1 = 2, M-l 
J = X 
K = 1 

TXEW(U.K) = T(1,J,K) + 2*FO*Rl*( T(EJ-l.K) - T(I.J,K) ) + 
& 2’ ; FO*R2*( T(I.J.K + 1) - T(I.J.K) ) + FO*( T(I+ 1.J.K) + 

& T(I-1.J,K) - 2*T(U.K) ) + 2*FO*BlO*(RI**0.5)*( TAMB - 
& T(I.J.K) ) 

170 CONTINUE 



4 PERIMETER NR. 10 



DO 175 1 = 2, M-l 
J = N 
K = P 

TXEW(I,J,K) = T(I,J,K) + 2*FO*RI*( T(I,J-1,K) - T(I,J,K) ) + 
& 2*FO*R2*( T(I,J,K-1) - T(I.J,K) ) + FO*( T(I + 1,J,K) + 

& T(I-I,J,K) - 2*T(U,K) ) + 2*FO*BIO*(Rl * !!! 0.5)*( TAMB - 
& T(I,J,K) ) 

175 CONTINUE 
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* PERIMETER NR. 11 



DO ISO 1 = 2, M-l 
J= 1 
K = P 

TNEW(I,J,K) = T(I,J,K) + FO*( T(I + 1,J,K) + T(I-1,J,K)- 
& 2*T(I,J,K) ) + 2*FO*Rl*(T(I,J + l,K)-T(I,J,K)) + 

& 2*FO*R2 <! ( T(I,J,K-1) - T(I,J,K) ) + 2* FO*Rl*DELY*FLUX/Kl 

180 CONTINUE 



* PERIMETER NR. 12 



DO 185 1 = 2. M-l 
J= 1 
K = 1 

TNEW(I.J,K) = T(I.J.K) + FO*( T(I+1.J,K) + T(I-1,J,K) - 
& 2*T(I J.K) ) + 2*FO*Rl*( T(I.J+ 1,K) - T(U.K) ) + 

& 2*FO*R2*( T(I J.K + 1 ) - T(I.J,K) ) + 2*FO*Rl*DELY*FLUX K1 

1S5 CONTINUE 



************************************************************************ 
**************** * * * * **************************************************** 

* UPDATING THE OLD NODE TEMPERATURES TO THE NEW 

* NODE TEMPERATURES. 

* 

DO 200 1 = 1, M 
DO 205 J= 1. N 
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DO 210 K= 1. P 



T(I.J.K)= TNEW(I,J,K) 

* 

210 CONTINUE 
205 CONTINUE 
200 CONTINUE 

4 

* 4 * 4 4 4 * 4 * * * 4 * * * 4 * * * * 4 * 4 4 4 4 * * * * * ❖ 4 4 * 4 * * * * 4 4 4 * * * * 4 4 4 4 * * * V- 4 * * * * * * * 4 * 4 4 * * * * 
4 $ $ 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 



* THIS PART OF PROGRAM DETERMINES IF SUBROUTINE VALID OR 

* SUBROUTINE VALID IS TO BE CALLED. IF THE USER DESIRES TO 

* VALIDATE THE PROGRAM IN CONJUNCTION WITH PROGRAM BRIAN 

* AND PROGRAM VALID. THEN THE SUBROUTINE OUTPUT WILL BE 

* BYPASSED. 

4 

IF(VAL .NE. 1) GO TO 250 
IF(TIME .EQ. DELT) CALL VALID 
I F( COUNT .NE. FREQ) GO TO 300 
CALL VALID 
GO TO 300 

A 

250 IF(TIME .EQ. DELT) CALL OUTPUT 
IF(COUNT .NE. FREQ) GO TO 300 
CALL OUTPUT 

* 

300 IF (TIME .LE. 3600) GO TO 1 

* NOTE: THE TIME PERIOD WAS ARBITRARILY TAKEN 

* AS 3600 SECONDS (1 HOUR). 

* 

PRINT* 
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PRINT* 



CALCULATING AND PRINTING CPU TIME 

TIME1 = CTIMEO 
CALL UMACH(2,NOUT) 

WRITE(NOUT,*) CPU TIME (SECONDS) = ', TIMEl-TIMEO 

STOP 

END 



* * « * * * * * * * * * * * * * * * * * * * * * * * * * $ * * * * * * « « * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
a * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * $ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 



* THIS SUBROUTINE WILL BE USED FOR PRINTING THE NODE 

* TEMPERATURES. 

SUBROUTINE OUTPUT 



DIMENSION T(25, 25.25). TNE\V(25, 25.25) 



REAL TAMB. II. RHO. ALFA. FO. FLUX, R. Rl. R2, K1 

REAL DELX, DELY, DELZ. DELT. TIME, TINIT. TNEW, BIO 

REAL TDIFE1. TDIFE2, TDIFE3 

REAL LX, LY. LZ 

REAL CTIME, TIMEO, TIME1 

INTEGER I, J, K, M, N, P, COUNT, FREQ, VAL, NOUT 

COMMON TNEW, TIME, DELT, M, N, P, TDIFEI, TDIFE2, TDIFE3, TINIT 
COMMON COUNT, FREQ 

WRITE(V) 'TIME = '.TIME 
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PRINT* 



350 DO 370 K= 1, P 

WRITER, *) 'TEMPERATURE DISTRIBUTION ON PLANE' 
WRITE(6,360) 'K-',K 
360 FORMAT(A3,12) 

PRINT- 

PRINT 380, ((TNE\V(I,J,K) f J = 1,N), I = 1,M) 

380 F0RMAT(1X,F6.1.1X,F6. 1,1X,F6.1,1X,F6.1,1X,F6.1,1X,F6.1,1X,F6.1, 
& 1X.F6. 1,1X.F6.1.1X,F6.1,1X,F6.1 

PRINT- 

PRINT- 

370 CONTINUE 

REINITIALIZING THE COUNTER TO ZERO 



COUNT = 0 



RETURN 

END 



* * * * * # * * * * * ft * * * * ft * * * * * * ft * * * * ft ft * * ft ft * * * * ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft 
ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft. ft ft ft ft ft ft ft ft ft. ft. ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft 



SUBROUTINE VALID 

DIMENSION T(25,25,25), TNEW( 25, 25.25) 

REAL TAMB, H, RHO, ALFA, FO, FLUX, R. Rl, R2, Kl 
REAL DELX. DELY, DELZ. DELT, TIME, TINTT, TNEW, BIO 
REAL TD1FE1, TDIFE2, TD1FE3 
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REAL LX, LY, LZ 

REAL CTIME. TIMEO, TIM El 

INTEGER I, J, K, M, N, P, COUNT, FREQ, VAL, NOUT 

COMMON TNEW, TIME, DELT, M, N. P, TD1FE1, TDIFE2, TDIFE3, TINTT 
COMMON COUNT, FREQ 

TDIFE1 = TNEW( 5,11, 5)- TINTT 
TDIFE2 = TNEW( 5, 5. 5)- TINTT 
TDIFE3 = TNEW( 5, 1, 5)- TINTT 

f* 

IF(TIME .NE. DE7 T) GO TO 450 
PRINT* 

PRINT* 

WRITE!*.*) 'TIME(SEC) TDIFE1 TDIFE2 TDIFE3' 

WRITE!*,' *) ' ' 

450 PRINT 460. TIME. TDIFE1. TD1FE2. TDIFE3 
460 FORMAT(lX.F$.3.3X.F7.3,3X.F7.3.3X,F7.3) 

•• REINITIALIZING THE COUNTER TO ZERO 



COUNT = 0 



RETURN 

END 



g g g g g g g * £ g g g g g g g g g g g g g g g g g g g g g g * g g a. * g ^ g * g g g g g g g * g g g * g g * g g * g * * * $ g * * g * g * <5 
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